A processing–modeling routine to use SNOTEL hourly data in snowpack dynamic models

Accepted Manuscript A processing-modeling routine to use SNOTEL hourly data in snowpack dynamic models Francesco Avanzi, Carlo De Michele, Antonio Ghezzi, Cristina Jommi, Monica Pepe PII: DOI: Reference:

S0309-1708(14)00128-6 http://dx.doi.org/10.1016/j.advwatres.2014.06.011 ADWR 2228

To appear in:

Advances in Water Resources

Please cite this article as: Avanzi, F., De Michele, C., Ghezzi, A., Jommi, C., Pepe, M., A processing-modeling routine to use SNOTEL hourly data in snowpack dynamic models, Advances in Water Resources (2014), doi: http:// dx.doi.org/10.1016/j.advwatres.2014.06.011

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

A processing-modeling routine to use SNOTEL hourly data in snowpack dynamic models Francesco Avanzia,∗, Carlo De Michelea , Antonio Ghezzia , Cristina Jommib , Monica Pepec a

Department of Civil and Environmental Engineering, Politecnico di Milano, Piazza Leonardo da Vinci, Milano, Italy b Department of Geoscience and Engineering, Delft University of Technology, Delft, the Netherlands c Institute for Electromagnetic Sensing of the Environment (IREA), Consiglio Nazionale delle Ricerche, Via Bassini 15, Milano, Italy

Abstract SNOTEL hourly and daily data are a strategic information about snowpack dynamics in western United States. Hourly data are highly noisy due to, e.g., non-physical temperature-based fluctuations of the signal or gauge under-catch. Noise may hinder, among other factors, the correct evaluation of precipitation events or the measurement of SW E, hence the reconstruction of accumulation and melt run-off timing. This makes hourly data practically useless without a denoising procedure. As this time resolution is widely used in hydrologic applications, here we test SNOTEL hourly data in modeling snowpack dynamics. A one-dimensional model of snow depth, snow water equivalent and bulk snow density has been adopted to this aim. We define an automated processing routine to denoise data-series of snow depth, snow water equivalent, bulk snow density, liquid and solid precipi∗

Corresponding author Email address: [email protected] (Francesco Avanzi)

Preprint submitted to Advances in Water Resources

June 28, 2014

tation. Special attention is paid into distinguishing the different types of precipitation, and processing snow depth data. Since sub-daily physical oscillations in snow depth data are difficult to be separated from instrument noise, a joint processing-modeling procedure has been designed. Forty SNOTEL sites throughout the western United States with multi-year data are considered for testing the procedure. The analysis shows that the model performance, expressed in terms of median values of the Nash-Sutcliffe coefficient, are higher than 0.8 for all the three variables, provided the first year of each dataset is used in the calibration phase. Keywords: Snowpack, mass dynamics, snow depth, SNOTEL, hourly data, data processing

1

1. Introduction

2

In the western United States, about 40% - 70% of precipitation falls in

3

the solid form [69], and around 70% of stream flow originates from snowpack

4

melting during spring and summer [27]. Therefore, the modeling of snow-

5

pack dynamics is strategic in assessing present and future scenarios of water

6

availability in the whole western North America [55, 11].

7

The SNOw TELemetry (SNOTEL) network has been collecting system-

8

atic data of snow water equivalent (SW E), total precipitation (P ), snow

9

depth (SD) and air temperature (T ) at an increasing number of weather

10

stations since the 60’s [69], permitting users to access a great amount of

11

data at daily and hourly resolutions. These also include data of bulk snow

12

density ρ, which can be derived by coupling SW E and SD data [32]. Sam-

13

pling duration (up to 40/50 years), the number of sites (roughly 800) and the 2

14

time resolution of measurements (daily or hourly) make SNOTEL network a

15

strategic source of data for hydrologic applications.

16

A typical SNOTEL site is equipped with a snow pillow, an ultrasonic

17

depth sensor, a thermistor and a rain gauge. Data from these instruments

18

are affected by many sources of uncertainty, such as missing values, out-of-

19

range readings, non-physical temperature-based fluctuations, wind-induced

20

vibrations of the snow depth sensor, or leaks from pillows and rain gauges

21

[15, 58, 30, 7, 2, 62].

22

SNOTEL daily data are subjected to a semi-automatic quality control

23

([69], SNOTEL staff personal communication, 2011), during and at the end

24

of each Water Year (W.Y.), to remove instruments noise. The high quality

25

of these data has made possible their use in a great number of studies, e.g.

26

[49, 69, 20, 55, 53, 27, 6, 22, 14, 23, 43, 51, 60]. Conversely, to our knowledge

27

no acknowledged processing routine is available for hourly data. At present,

28

these are made available to the public in the form of raw data-series, which

29

strongly limits their use in hydrologic applications, since it hinders: 1) the

30

interpretation of cumulative precipitation data, i.e. the reconstruction of

31

solid and liquid events, which play a key role as input variables in hydrolog-

32

ical models (see [7] as an example), 2) the evaluation of snow depth and its

33

duration on the ground, useful for remote sensing applications, or 3) the cor-

34

rect registering of SW E dynamics, hence the reconstruction of accumulation

35

and melt run-off timing. As a result, very few contributions using SNOTEL

36

hourly data are available in the literature, such as [7, 3, 63, 16].

37

Nonetheless, high-resolution data are essential in many hydrologic ap-

38

plications [26, 17], and especially in modeling the dynamics of snowpacks

3

39

[72, 8, 9, 16, 39]. In particular, the hourly resolution allows 1) capturing the

40

day-night cycles of melting and refreezing [59, 47, 73, 5], 2) displaying the

41

oscillations of air temperature and their relation with the snowpack outflow

42

[67, 75] and 3) evaluating the dynamics of precipitation, the principal mass

43

input of the seasonal snowpack [44, 4]. As a matter of fact, modern snowpack

44

models operate on a sub-daily, or even sub-hourly, time step (see for exam-

45

ple SNTHERM, [33], CROCUS [12, 74, 13], UEB [72, 45, 46], SNOWPACK

46

[8, 41, 42, 75], ALPINE3D [9, 10] or the model proposed by [36]), and many

47

automatic weather stations all around the world are nowadays equipped to

48

register data at the hourly resolution. Examples are the Col de Porte site in

49

Grenoble, France [54], the Weissfluhjoch site in Davos, Switzerland [8, 66],

50

or the JMA network in Japan [76]. The coupling of snowpack models and

51

high-resolution data is useful in a wide set of applications, such as avalanche

52

prediction, remote sensing imaging interpretation and runoff dynamics as-

53

sessment [75]. It follows that an acknowledged data processing routine for

54

SNOTEL hourly data could contribute to better exploit in the future this

55

huge body of data for snowpack modeling applications.

56

Here, we propose a simple and automated data pre-processing routine for

57

SNOTEL hourly data of SW E, T , SD, P and ρ. Because of the difficulty

58

in discerning between sub-daily oscillations of the variables (due to processes

59

like snow and rain events, melting and settling) and noise, the routine has

60

been coupled to a model [16] (hereinafter, HyS model), to assist in the de-

61

noising of the signal. We refer to this approach as the processing-modeling

62

routine. We test this routine extensively using 40 SNOTEL sites.

4

63

64

2. Data An exhaustive presentation of SNOTEL network, including its climatol-

65

ogy, can be found in [69]. Data are available at the URL

66

http://www.wcc.nrcs.usda.gov/snow/.

67

2.1. The SNOTEL network: site instrumentation

68

At each SNOTEL site, all measurements are automatically logged. A

69

snow pillow collects the current snow water equivalent, SW E, by measuring

70

through a pressure transducer the weight of the snow cover, which overlies

71

a plate located on the ground surface. An ultrasonic depth sensor provides

72

the measurement of snow depth on the ground, SD. The sensor is located

73

above the ground at a suitable height, and measures the time of flight needed

74

by a 50 kHz signal to travel to and from the snow top surface below. By

75

evaluating the speed of sound in air as a function of air temperature, the

76

sensor is able to detect its distance from the snow surface, and therefore to

77

derive SD. Total precipitation is held inside a heated rain gauge, a cylindrical

78

vessel, from 2.5 m to 5 m tall, with a 0.3 m orifice. It does not distinguish

79

between solid and liquid precipitation. Air temperature is measured with a

80

thermistor. The resolution of the instrumentation is 0.1 inch (2.54 mm) for

81

SW E and precipitation, 1 inch (2.54 cm) for snow depth, and 0.1◦ C for air

82

temperature.

83

2.2. The SNOTEL network: data uncertainty

84

85

The main sources of uncertainty in the measurement at snow pillows include [15, 19]:

5

86

• Transmission errors, including missing values (flagged with a -99.9

87

value) and random out-of-range readings, i.e. extremely high or nega-

88

tive values;

89

• Temperature-based fluctuations of the signal (sometimes known as flut-

90

ter in SW E), caused by changes in air temperature, affecting the den-

91

sity of the antifreeze solution inside the pressure transducer, hence its

92

stress state [50, 70]. Also, similar effects are ascribable to the possible

93

presence of air bubbles in pressure transducer lines [15]. The absolute

94

value of the fluctuations is usually low, comparable to the instrument

95

resolution;

96

• Slope effects, as creep [61, 35];

97

• Over/under measurements due to the thermal and mechanical differ-

98

ences between the pillow, the snow and the surrounding soil [29, 30,

99

31, 34]. These include differential melting, snow bridging and ice lenses

100

intrusions;

101

• Instrumental malfunctioning, including leaks [15].

102

Precipitation data are affected by:

103

• Transmission errors, including missing data (flagged with a -99.9 value)

104

and extremely high or negative readings (as for the snow pillow);

105

• Instrumental under-catch, caused by wind turbulence and plugging [7];

106

• Instrument malfunction, mainly due to evaporation and leaks.

6

107

108

• Temperature-based fluctuations, as for the snow pillow, since the measurement is taken through a pressure transducer

109

They all prevent accurate registration of many events, returning, in gen-

110

eral, an underestimated recording of precipitation [40, 24, 18]. This is par-

111

ticularly evident when solid events occur, with a catch deficiency up to 50%

112

- 70% of the total snow precipitation, depending on wind velocity [40].

113

For temperature data, transmission errors represent the only evident

114

source of uncertainty. Thermistors also can be affected by unphysical heating

115

during the day, or cooling during the night, if not installed properly [58].

116

117

Few studies are available in the literature about the uncertainties of ultrasonic depth sensors readings [25, 71, 58, 62, 2]. These include:

118

• Transmission errors;

119

• Interferences within the Field-of-View (FOV) of the instrument due to

120

the transit of snowflakes and droplets, or caused by grass at the border

121

of the measurement area, when SD is small;

122

• Wind-induced sensor vibrations [62];

123

• Erroneous sub-daily data fluctuations [58], also known as flutter in SD.

124

These are partly due to temperature oscillations, which cause a varia-

125

tion in the speed of sound (SOS) of the sensor signal, and to the effects

126

of temperature on the electronics of the sensor [58]. The instrument

127

measures T to correct the value of SOS, but fails in completely fix-

128

ing hourly erroneous oscillations. This kind of data uncertainty causes

129

repeated oscillations of recorded snow depth, around the actual value. 7

130

The absolute value of these oscillations is comparable to the instrument

131

resolution, i.e. 2.54 cm.

132

The overlap of all these sources of uncertainty makes raw data noisy and

133

their use in the study of snowpack dynamics at the hourly time scale difficult.

134

This could explain the rare use of SNOTEL hourly data in snow hydrology

135

literature.

136

2.3. The considered SNOTEL sites

137

Here, we have considered 40 out of ∼ 800 SNOTEL sites (i.e., ∼5%),

138

chosen according to the following criteria: 1) each State of western U.S. has

139

to be represented by at least one site, 2) the elevation of the sites has to be

140

representative of the elevation range of SNOTEL network, 3) datasets have

141

to be complete. Table 1 reports the sites, and the water years considered.

142

In Figure 1(a) and Figure 1(b), sites are depicted on Landsat TM images.

143

Figure 2 gives the histogram of sites elevation. Half of the sites are placed

144

between 2000 and 2800 m.a.s.l., due to the altitude range of the gauged

145

area, but lower (< 2000 and > 1200 m.a.s.l.) and higher (> 2800 m.a.s.l.)

146

elevations are also represented. The two sites at an altitude lower than 1200

147

m a.s.l. are situated in Alaska. As a whole, we have considered 268 years,

148

with an average of 6.7 years per station.

149

3. Methods

150

151

In this section, we present the proposed processing routine, the model used, and their coupling in a processing-modeling routine.

8

152

3.1. The data processing routine

153

In this section, we propose a processing routine to denoise hourly data of

154

air temperature, precipitation, snow water equivalent and snow depth. All

155

the processing steps are based on an acceptance-rejection approach, i.e. a

156

threshold-based definition of anomalous values. Here, we propose a rejection

157

criterion, implying the replacing of the value at a certain hour with the one

158

at the previous hour. Replacing the missing value is necessary to couple the

159

data with snowpack models. This assumption works well for SW E and SD

160

variables, usually characterized by a slow variation rate, and less well for T

161

and P data, having a quicker variation rate. Figure 3 shows all the tasks of

162

the data processing. The processing routine is a stand-alone procedure, i.e.

163

it corrects the data from each site using only the information of that site.

164

SW E and P data are processed in mm, SD in m, and T in

◦

C. Hence,

165

a preliminary step of data conversion from US customary units is needed.

166

3.1.1. SW E and T data processing

167

The processing routine firstly considers SW E and T data, and removes

168

transmission errors. The other sources of uncertainty are not processed here

169

either because of their small magnitude (e.g. flutter in SW E), or because

170

no proxy indication which can assist in this operation is available (as for the

171

over/under-measurements of snow pillows, or non-physical heating/cooling

172

of the thermistor). Creep effects are not considered since SNOTEL sites are

173

usually installed on an horizontal surface.

174

Negative values and data which exceed a fixed threshold, SW EM AX (usu-

175

ally 150 inches, i.e. 3800 mm, but it can be site-dependent), are rejected. To

176

our knowledge, this step is similar to the imposition of lower and upper limit 9

177

profiles usually operated by SNOTEL staff on daily data. In addition to this,

178

the procedure sets at zero any SW E value recorded during summer. The

179

snow season is defined, at each site, through the first (MM IN ) and the last

180

(MM AX ) month. Typically, MM IN = 10 (October) and MM AX = 6 (June)

181

are assumed, but they can be site-dependent. The definition of summer by

182

means of fixed temporal parameters is preferable to other strategies, such

183

as zeroing SW E when SD is null, and/or when air temperature is stably

184

positive for a very long period, say some weeks, since it does not depend on

185

other measured variables, which can be subjected to noise too.

186

As for T , the value at a given hour is rejected if its absolute value is greater

187

than a specific threshold T M AX , usually assumed equal to 40◦ C. T M AX can

188

be calibrated at each site.

189

3.1.2. SD data processing

190

The procedure involves a set of basic tasks used to correct transmission

191

errors. These are quite similar to those proposed by [69] and [52], and those

192

sometimes adopted to our knowledge by SNOTEL staff in providing daily

193

data. They are 1) rejection of any SD value, which is negative, or greater,

194

than a specific threshold (HM AX , site-specific); 2) zeroing any SD value

195

reported during summer. The SD data, at the end of these steps, are stored

196

as SDRAW .

197

Secondly, the procedure removes the erroneous sub-daily fluctuations in

198

SD data. This operation is necessary because 1) fluctuations in SD data are

199

much larger than the ones in SW E data; 2) refined SD data can be useful to

200

determine solid precipitation amounts [65]. The routine computes the hourly

201

differences (J) of SD, and rejects each increment (/decrement) equal to the 10

202

absolute value of a decrement (/increment) in the previous, or successive,

203

time step, to remove trivial readings errors. In addition, the routine rejects

204

every SD increment (/decrement) larger (in modulus) than a fixed threshold

205

DHM AX (usually fixed at 0.6 m). This step is similar to the imposition of

206

lower and upper limit for daily increments(/decrements) sometimes operated

207

by SNOTEL staff on daily data (personal communication, 2011). Then, a re-

208

fined step is established to reduce SD oscillations due to temperature-based

209

fluctuations. We assume as correct any SD variation which is associated

210

with stable air temperature conditions, and as erroneous any SD variation

211

accompanied by abrupt temperature variations. This is in agreement with

212

experimental observations [58]. On the contrary, no acknowledged procedure

213

is known to fix hourly oscillations. To this aim, the absolute value of the

214

difference between air temperature and the average air temperature of the

215

previous 24 hours, i.e. ∆T24 (t), is computed for each hour. Any SD value

216

associated with a ∆T24 (t) greater than a threshold (∆TM AX ) is rejected.

217

Thus, SD does not vary until a small fluctuation in air temperature is reg-

218

istered. In other words, ∆TM AX is an upper threshold of temperature stable

219

conditions.

220

As a third step, J is re-computed, and the positive increments are gath-

221

ered to form the input array of solid events (s). To avoid any misconstruction

222

of events due to interferences in the field of view, which can occur when SD

223

is small, every value of s registered during the last part of the melting season

224

(fixed by the definition of the last month of possible snow events, MLAST ,

225

usually April), is removed.

226

At the end, the cumulative solid precipitation curve, H, is calculated, year

11

227

by year, in mm of water equivalent. To this aim, density of fresh snow, ρN S ,

228

is calculated using the formula proposed by [1]. Accordingly, ρN S [kg/m3 ]

229

= 50+1.7 (T + 15)1.5 , if T is greater than -15 ◦ C, and 50 kg/m3 otherwise. In

230

multi-year processing, the cumulative curve of solid precipitation is zeroed

231

every 1st October. In this way, the ultrasonic depth sensor is used as the

232

unique predictor of solid events, overcoming the well-known difficulties of

233

rain gauges in registering these events, because of plugging, and wind effects.

234

3.1.3. P data processing and solid/liquid partitioning

235

Filtering of P data is necessary to reconstruct liquid events. Similar to SD

236

processing, preliminary operations include 1) the rejection of any negative,

237

or erroneously high, value (through the definition of the threshold PM AX );

238

2) the rejection of any P increment (/decrement) larger than a threshold,

239

DPM AX . These operations fix transmission errors. Then, 3) any decrement

240

is rejected. These preliminary operations are quite similar to those usually

241

operated by SNOTEL staff on daily data of cumulative precipitation.

242

To determine liquid events, the routine calculates the daily variations of

243

P , and compares them with the contemporaneous variations in the cumula-

244

tive solid precipitation curve, H. Every positive difference between P and

245

H is considered as a liquid event, and reported in the array, L, of liquid

246

events at the end of the day. Solid and total events are derived from two

247

independent instruments, and some incongruities can occur. In winter, the

248

daily increment of solid precipitation is often greater than the correspond-

249

ing increment of the total cumulative precipitation, because of the afore-

250

mentioned instrumental under-catch of rain gauge. In this case, no further

251

liquid contribution is added. The procedure considers the rain gauge as a 12

252

reliable instrument for the recording of liquid events, as described by [40].

253

This mass-conservation approach to solid-liquid partitioning is preferred to

254

proxy-based methods (e.g. temperature-based) because at SNOTEL gauges

255

different sensors are available to monitor precipitation, and because of the

256

site-specificity of proxy-based methods [37]. Each hourly positive value of

257

liquid precipitation is indicated with p.

258

At the end of the procedure, a cumulative precipitation curve, PM OD ,

259

is calculated by summing H and L contributions. We use the daily time

260

resolution as control window of precipitation to avoid local incongruities in

261

the data. Among them, for example, delay in registration of solid events may

262

occur due to the sticking on vessels walls [7]. The definition of a routine at

263

the hourly scale is an important future development in this framework.

264

3.1.4. Snow density evaluation

265

Bulk snow density (ρ) values are obtained by coupling edited values of

266

SW E and SD, only when SD > 0. To avoid non-physical values, which can

267

be computed when SW E and/or SD values are low, ρ values larger than

268

1000 kg/m3 are rejected. This because 1000 kg/m3 is the maximum density

269

of the ice-air-liquid water mixture.

270

At the end of the data processing, edited arrays of SD, SW E, T , P and

271

ρ are available, together with arrays of cumulative precipitation in the two

272

forms, solid and liquid.

273

The proposed approach to data processing requires the definition of some

274

convenient

thresholds:

275

SW EM AX ,

276

DPM AX , PM AX , ∆TM AX and MLAST .

MM IN ,

T M AX ,

MM AX ,

13

HM AX ,

DHM AX ,

277

Some of these (SW EM AX , MM IN , MM AX , T M AX , HM AX , DHM AX , DPM AX ,

278

PM AX ) can be easily calibrated by visual inspection of data. These thresh-

279

olds have a climatological base and can be estimated for each site once, and

280

should not change throughout the years.

281

The evaluation of ∆TM AX and MLAST cannot be performed by visual in-

282

spection, because flutter masks the real increases and decreases of SD among

283

many non-physical fluctuations, and because it is impossible to define a pri-

284

ori 1) the maximum value of ∆T24 (t) that does not affect erroneously SD

285

data, and/or 2) when FOV interferences begin affecting SD measurements.

286

Thus, their calibration is operated by coupling the processing routine to a

287

model, which uses the processed arrays as inputs, as in [16]. In this way, er-

288

roneous input events can be detected as physical incongruities between data

289

and model.

290

3.2. The processing-modeling approach to parameters calibration

291

3.2.1. The HyS model

292

The HyS model [16] has a very simple formulation so that it can be

293

used in conditions of data scarcity, and with low computational efforts. It

294

aims at predicting SW E dynamics at a point through the coupled modeling

295

of bulk snow density and snow depth, by considering the snowpack as an

296

equivalent one-layer two-constituents (one dry and one wet) mixture. The

297

dry constituent includes ice and air, while the wet one includes liquid water

298

only. The model uses the depth of the ice structure, hS , its dry density ρD ,

299

and the depth of the liquid constituent, hW as state variables. The water

300

density ρW is fixed at 1000 kg/m3 .

301

The dynamics of snowpack, at the site scale, are computed by imposing 14

302

mass conservation of the two constituents, momentum conservation of the

303

dry phase (which is the constituent subject to deformations) and assuming

304

a viscous constitutive equation. The model variables and parameters are

305

reported in Table 2. Table 3 summarizes the governing equations, which

306

are integrated with a time step of one hour to describe snowpack dynamics.

307

The considered mass fluxes are solid and liquid precipitation, melting and

308

water outflow. All the other contributions (such as the fluxes due to wind,

309

sublimation, evaporation, refreezing and soil infiltration) are neglected. Only

310

two model parameters must be calibrated, namely the degree-hour coefficient

311

a and the runoff parameter c. The input variables are s, p and T .

312

This model does not account for melting-refreezing dynamics, sub-daily

313

variations in the energy balance and the effect of wind on snow accumula-

314

tion/depletion. Nonetheless, it is simple and physically-based, hence useful

315

as a first tool to test the feasibility of using hourly SNOTEL data in snow-

316

pack modeling. This contribution is a first step in this direction, while fu-

317

ture attempts should verify the quality of hourly data when using them in

318

physically-based snowpack models with detailed parametrization.

319

3.2.2. The processing-modeling routine

320

The processing routine, described in 3.1, needs some parameters to be

321

calibrated. Among these, ∆TM AX and MLAST cannot be calibrated sepa-

322

rately from a model. On the other hand, the considered model needs the

323

calibration of two parameters, a and c. As data processing and snowpack

324

modeling are coupled, ∆TM AX , MLAST , a and c are calibrated jointly, by

325

minimizing the errors in a simulation.

326

For each site, using the entire dataset, we initially estimate SW EM AX , 15

327

MM IN , MM AX , T M AX , HM AX , DHM AX , DPM AX , PM AX through visual in-

328

spection. Then, we evaluate MLAST , eliminating local overestimations of

329

modeled snow depth at the end of the melting season. Notice that this last

330

operation is marginal in reproducing the whole season.

331

Secondly, the least-square method is used between processed SD data and

332

modeled h (control depth, see Table 2), to evaluate the optimum parameters

333

set (∆TM AX , a), for the considered site, and the year of calibration. Hence,

334

the minimum of the error surface, quantified using square differences between

335

observations and model, in the ∆TM AX - a plane is searched.

336

Thirdly, the least-square method is applied to observed and modeled val-

337

ues of snow density, ρ, to evaluate c. Although c can potentially vary from

338

0 to +∞, its range of variation is restricted to (0 ÷ 1], because of numerical

339

instabilities [16].

340

3.2.3. Evaluation of the procedure

341

The evaluation of processing-modeling performances is conducted as fol-

342

lows: 1) we first assess the effect of the processing routine on raw data. The

343

anticipated outcome is a processing routine which removes noise, without sig-

344

nificantly affecting the measurements; 2) we secondly operate a calibration of

345

the whole set of parameters, and evaluate procedure performances by means

346

of Nash-Sutcliffe coefficients NSE [56, 48, 64], comparing modeled results and

347

processed data-series. High values of average NSE (say, greater than 0.6),

348

irrespective to the calibration year, would denote a successful use of hourly

349

SNOTEL data in snowpack modeling; 3) thirdly, the variability of parame-

350

ters, nivometric coefficients and catch-ratio estimations is commented.

351

As for the calibration, we have initially considered a first-year calibration 16

352

procedure, i.e. we have applied the processing-modeling routine to the first

353

year available of each dataset, and evaluated the performances of the rou-

354

tine by using the Nash-Sutcliffe coefficient NSE, calculated for each of the

355

variables SW E, SD and ρ by comparing the processed data with the mod-

356

eled ones. The coefficient has been calculated both on the calibration year

357

(indicated as NSESD , NSEρ and NSESWE respectively for SD, ρ and SW E),

358

and on the entire dataset as average of the year values (indicated as NSESD ,

359

NSEρ and NSESWE ), excluding summer periods.

360

As a second approach, we have investigated how routine performances and

361

parameters values vary, when changing the calibration year. We have used

362

one at the time all the available years in the calibration phase, and we have

363

evaluated the performance, year by year, on the whole dataset, by means of

364

the NSE coefficients. We indicate this as one-year calibration procedure.

365

4. Results and discussion

366

During the calculations, 13 years of the dataset have been discarded (∼

367

4.8%), because of instrument errors. Thus, the actual dataset has a total

368

duration of 255 water years, with an average of 6.3 years/site.

369

4.1. Raw data processing

370

Here, we consider SNOTEL data at site S26 (randomly chosen) as a case

371

study to illustrate the effects of the processing routine on raw data (point 1,

372

Section 3.2.3).

373

In Figure 4, we report the comparison between raw (green) and processed

374

(black) SNOTEL hourly data during the Water Year 2011. In particular,

375

panel (a) shows the comparison in terms of snow depth, while panel (b) 17

376

reports raw and processed cumulative precipitation data. Panel (b) gives

377

also processed values of cumulative solid and liquid precipitation. In panel

378

(c), the comparison in terms of SW E is reported, while in panel (d) a zoom

379

of panel (a) is given for sake of clarity. Finally, we report in panel (e) raw and

380

processed hourly data of precipitation (in both the forms). Data have been

381

obtained by selecting hourly positive differences of the raw and processed

382

cumulative precipitation data-series, respectively. Raw data in Figure 4

383

have been already filtered for negative and anomalously high readings.

384

Summer erroneous readings are very evident in panel (a), as also flutter

385

oscillations. These oscillations would lead to overestimations of modeled

386

snowpack depth and SW E, if not processed, since every positive increment

387

would be considered as an event. A specific focus on flutter oscillations is

388

reported in panel (d). It shows the effect of flutter at the sub-daily resolution,

389

and the smoothing operated by ∆TM AX . On the contrary, the proposed

390

processing routine has a marginal effect on SW E data (panel c): processed

391

data differ from the raw ones during summer, but the two data-series coincide

392

during the rest of the W.Y.

393

Regarding P data, panel (b) shows 1) the correction of catch deficiency,

394

visible as the difference between the raw P dataset and the processed one,

395

and 2) the partition of events. Panel (e) shows also the great amount of

396

fictitious precipitation events one would erroneously infer by simply assuming

397

that each positive difference in raw data is an event. These are due to the

398

occurrence of oscillations in raw data. The sum of these positive differences

399

on W.Y. 2009 is equal to ∼ 11500 mm, i.e. 4.6 times greater than the

400

processed total precipitation at the end of the W.Y. (see panel (b)).

18

401

The processing-modeling routine does not alter the general shape of SD

402

data-series. This is shown in Figure 4, where the smoothing of flutter oscilla-

403

tions is evident, but preserving the global evolution of SD at the same time.

404

Similarly to a Nash-Sutcliffe coefficient, NSER has been calculated between

405

processed SD and SDRAW for every calibration year, to quantify the impact

406

of flutter and reading errors corrections on the reconstruction of the raw SD

407

data-series. This coefficient is greater than 0.985 in 239 out of 255 years

408

(∼ 94%) considered during the calibration phase. Thus, processed and raw

409

data-series are very similar to each other, apart from second-order differences

410

(i.e., mainly flutter).

411

Generally speaking, the routine succeeds in maintaining the shape of data.

412

This has been visually verified on each of the 255 W.Y. considered.

413

4.2. Evaluation of routine performances

414

Here, we focus on the evaluation of processing-modeling routine perfor-

415

mances when simulating SW E, snow depth and bulk snow density (point 2,

416

Section 3.2.3).

417

Figure 5 reports an example of processing-modeling routine results at site

418

S33 (randomly chosen) for the period 2006-2011, obtained using the first-year

419

calibration procedure. In panel (a), modeled time series of h (red) are plotted

420

vs. processed SD data (black). Raw data are not reported. This comparison

421

shows the agreement between snow depth data and model, during both the

422

accumulation and the melting seasons, even using just one year (the first one)

423

in the calibration phase, and five years in the validation. In addition, hW

424

time series (blue) are reported. This variable increases during the melting

425

season, as expected. 19

426

In panel (b), processed bulk snow density (black), modeled ρ (red), and

427

modeled ρD (blue) time series are displayed. The modeled variables are close

428

to each other throughout the accumulation season, while during the melting

429

season ρ tends to increase because of liquid water saturation. This is in

430

agreement with data.

431

In panel (c), processed (black) and modeled (red) SW E are shown. Their

432

good agreement, during the first year of calibration and the following five

433

years of validation, illustrates the performances of the routine.

434

We report in Figure 6 a comparison at site S26 (W.Y. 2009) between

435

processing-modeling results (panel (b), (d) and (f) for snow depth, snow

436

density and SW E, respectively) and the modeled results one would obtain

437

by using directly raw input data (raw-modeled results, panel (a), (c) and

438

(e) for snow depth, snow density and SW E, respectively), and the set of

439

parameters obtained from the processing-modeling routine. In the panels,

440

green lines represent raw data. The comparison between panel (a) and (b)

441

shows the great overestimation of snow depth caused by summer and flutter-

442

induced fluctuations. These are evident considering 1) the erroneous increase

443

in raw-modeled data at the beginning of the year, result of summer erroneous

444

readings, 2) the constant increase in h during the period between 31/12/2008

445

and 1/4/2009, when raw-modeled data return a total increase equal to ∼ 600

446

cm, while raw data measured a global increase of snow depth of ∼ 70 cm

447

(i.e., ∼ 11 %), and 3) the erroneous persistence of a positive snow depth at

448

the end of the year (∼ 5 m). In the same panel, it is reported in purple the

449

modeled snow depth one would obtain using data processed using just the

450

visual thresholds (i.e., without considering ∆TM AX and MLAST ). Much of

20

451

the overestimation during the period of existence of the snowpack is due to

452

erroneous sub-daily data fluctuations. The same considerations apply also

453

to SW E (panel (e) vs panel (f)), which is strongly overestimated using raw

454

input data, and even using visual thresholds only (purple line). Overesti-

455

mation of bulk snow density is moderate (panel (c) vs panel (d)), regardless

456

of the degree of processing used. This example demonstrates that data pro-

457

cessing is mandatory to use raw SNOTEL hourly data in a model, and that

458

the processing-modeling routine succeeds in removing large amounts of noise,

459

mainly induced by erroneous sub-daily data fluctuations.

460

4.2.1. The first-year calibration

461

Figure 7 summarizes the results of the first-year calibration procedure.

462

NSE values have been calculated both on the calibration year (left), and on

463

the entire dataset (right). The negative values of the Nash-Sutcliffe coeffi-

464

cients are set to 0. In the figure, crosses represent outliers.

465

For the calibration year, NSE values are very high, especially for SD.

466

This because two parameters (∆TM AX and a) are involved in the processing-

467

modeling of SD, while just one parameter (c) is used to model bulk snow

468

density dynamics. The median values of NSE are ∼ 0.945 for SD and ∼

469

0.90 for SW E and ρ. The maximum values of NSE are ∼ 0.99 for all the

470

variables, while the minimum values are ∼ 0.85 for SD, ∼ 0.62 for SW E

471

and ∼ 0.7 for ρ.

472

NSE values are high, but have a wider variability. The minimum NSESD

473

coefficient is ∼ 0.6, and three outliers are present. The median value is equal

474

to 0.85, while a maximum value of 0.95 is reported. As for NSESWE , the me-

475

dian/maximum/minimum values are equal to 0.78/0.96/0.45, respectively, 21

476

while four outliers are reported. On the contrary, the median/maximum/minimum

477

values of NSEρ are equal to 0.85/0.95/0.63, respectively.

478

SW E performances are affected by the performances of both SD and ρ

479

simulations. In fact, they present the widest range of NSE. On the contrary,

480

NSEρ preserves the excellent performances of the calibration period. As a

481

whole, the processing-modeling routine succeeds in returning values of NSE ≥

482

0.6 (see Section 3.2.3), even using just one year (i.e., the first one available)

483

in the calibration phase.

484

4.2.2. The one-year calibration

485

The results of the one-year calibration are reported in Figures 8. Panels

486

(a), (b) and (c) show respectively the box-plots of NSESD , NSEρ and NSESWE

487

for the considered sites. Negative values of the Nash-Sutcliffe coefficient are

488

set to 0 for graphical purposes, and outliers reported as crosses.

489

The minimum values of the Nash-Sutcliffe coefficients are greater than

490

0.6 at 27 out of 40 sites (∼ 67%) for NSESD , 34 sites (∼ 85%) for NSEρ , and

491

19 sites (∼ 47.5%) for NSESWE . Furthermore, minimum values of NSE are

492

greater than 0.8 at 13 sites (∼ 32.5%) for SD, 23 sites (∼ 57.5%) for ρ, and

493

6 sites (∼ 15%) for SW E. Hence, the routine performances are successful

494

throughout the western United States, even using a one-year procedure for

495

the calibration.

496

In particular, the bulk snow density is the variable with the highest values

497

of NSE. This observation is in agreement with the findings reported in the

498

previous section. The width of the boxplots in panel (a), and the presence of

499

outliers characterize the variability of processing-modeling performances on

500

SD. In some cases, such as sites S6, S9, S13, S18 and S25, the calibration 22

501

year is crucial in obtaining good performances. Nonetheless, since at all

502

sites the median value of NSESD is ≥ 0.5, low performances seem to be

503

isolated, probably due to instrumental reasons. Results are more variable

504

when considering NSESWE . At 16 sites, at least one value of NSESWE is ≤ 0.

505

As for S6, S9, S13, S18, S23 and S25, the width of NSESWE box plot is mainly

506

ascribable to SD performances, while at S3, S15, and S40, the low values of

507

NSESWE are due to low performances in bulk snow density modeling.

508

4.3. Parameters variability

509

Here, we comment on the variability of parameters’ values, nivometric

510

coefficients and rain gauges catch deficiencies (point 3, Section 3.2.3).

511

4.3.1. ∆ TM AX and a

512

Panel (a) of Figure 9 reports an example of the square error surface

513

in ∆TM AX - a plane, calculated for S4 site and the water year 2004 using

514

SD processed data and the modeled ones. Panel (b) of Figure 9 shows the

515

contour lines of the same error surface. The surface shows that big errors

516

are associated to two possible combinations: 1) low values of ∆TM AX and

517

high values of a, and 2) high values of ∆TM AX and low values of a. On the

518

contrary, low square error values are reported within two areas, identified by

519

1) low values of ∆TM AX and low values of a, and 2) high values of ∆TM AX

520

and high values of a. These two areas are connected by a broad region of

521

error equivalence, sub-parallel to the main diagonal of panel (b). In this case,

522

the optimum couple is ∆TM AX = 2.2◦ C and a = 0.00026 m/h/◦ C.

523

The shape of the error surface is explained by considering that ∆TM AX

524

determines directly the number of snow events, by affecting the variability of 23

525

SD and, therefore, the presence of any s. A low value of ∆TM AX causes the

526

rejection of a relevant part of SD oscillations resulting in a flat signal in the

527

limit case ∆TM AX → 0. Conversely, high values of the same parameter (say,

528

∆TM AX ≥ 3◦ C) tend to retain the most part of SD oscillations resulting in

529

a high variable signal. A high value of a implies a rapid melting of the ice

530

structure, and therefore a rapid decreasing in hS , while on the contrary a

531

low value of a would imply an opposite behavior, given identical atmospheric

532

conditions. This confirms the coupling of the two parameters, and the need

533

for a joint optimization.

534

Figure 10(a) reports the box-plot of a parameter, for the 40 sites, using the

535

one-year calibration. The external (i.e., intersite) variability of the parameter

536

is wider than the internal (i.e., intrasite) variability, with very few outliers.

537

The a parameter is therefore quite site-specific, in agreement with previous

538

analysis reported in [28].

539

Figure 11(a) shows the variability of median values of a with altitude.

540

The parameter tends to slightly increase with site elevation. Again, this is

541

in agreement with the analysis by [28]. Nonetheless, the elevation fails in

542

completely predicting the a variability, probably because of its dependence

543

from other factors (e.g. the season).

544

Figure 10(b) reports the boxplot of ∆TM AX for the 40 sites. In 27

545

cases out of 40, the median value of this parameter is smaller or equal than

546

2 ◦ C. This means, conversely, that at most sites a ∆T24 (t) greater than

547

2 ◦ C will lead to a SD erroneous oscillation. This confirm the noisiness

548

of SD data. Many sites - roughly 14 out of 40 - have a variation range

549

which is roughly equal, or greater than, 4 ◦ C. In many of these, NSESD

24

550

(and NSESWE ) has 1) a wide variation range; 2) low values (Figure 8).

551

This confirm the importance of data pre-processing and the sensitivity of the

552

routine performance to ∆TM AX calibration.

553

4.3.2. Runoff parameter c variability

554

Figure 10(c) reports the boxplot of c for the 40 sites. At 25 sites out of 40,

555

the calibration values vary over roughly 60% of the possible variation range

556

(0 − 1]. This great variability has no elevation gradient, as shown in Figure

557

11(b), and it can be explained by considering that c variability is affected by

558

1) hW dynamics; 2) ρD dynamics, which have been taken from literature and

559

not calibrated; 3) the complicated water pattern in snow [68]. Nonetheless,

560

the processing-modeling routine is not particularly sensitive to c variability,

561

since it hardly influences the values of NSEρ and NSESWE (Figure 8).

562

Figures 10(c) and 11(b) show a concentration of c values equal to 1.

563

This is due to the fact that the parameter drives the hydraulic dynamics

564

of snowpack melting, and rules the distance between ρD and ρ. Thus, an

565

overestimation of ρD (not calibrated) leads necessarily to an overestimation

566

of c and, in some cases, causes the observed values of c = 1.

567

4.3.3. Solid/liquid partitioning

568

Figure 12(a) reports the nivometric coefficient of each year as a function

569

of the site elevation. It was calculated using the optimum couple (∆TM AX ,

570

a) obtained by the calibration routine for that year. The median value of all

571

the nivometric coefficients is equal to 0.62, the minimum one to 0.23, and

572

the maximum one to 0.93. The first and the third quartiles are 0.53 and 0.7,

573

respectively. These results are in agreement with the values given by [69] for 25

574

the western United States.

575

In Figure 12(b), the nivometric coefficient of each calibration year has

576

been plotted as a function of the correspondent ratio PM OD /P between the

577

modeled and observed values of cumulative precipitation at the end of each

578

calibration year. Again, it has been calculated using the optimum couple

579

(∆TM AX , a) obtained by the calibration routine for that year. Results show

580

that the higher the nivometric coefficient, the higher is this ratio. The median

581

value of PM OD /P is ∼1.5. The first and third quartiles are, respectively,

582

∼1.4 and ∼1.6. Conversely, median and quartiles values of P/PM OD (gauge

583

undercatch, [21]) are respectively ∼0.67, ∼0.71 and ∼0.61. These values

584

are in agreement with the indications given by [40], supporting the results

585

of the processing-modeling routine both in 1) reconstructing the physical

586

cumulative precipitation curve, and in 2) partitioning precipitation events.

587

5. Conclusions and Outlook

588

We have discussed the use of SNOTEL hourly data in one-dimensional

589

snowpack modeling. To this aim, a simple model of the temporal dynam-

590

ics of snow water equivalent, snow depth and bulk snow density (Hys model,

591

[16]) has been used, and a processing routine has been introduced, to denoise

592

data of snow depth, snow water equivalent, air temperature and precipita-

593

tion, partitioned in solid and liquid events. This routine is based on simple

594

climatological criteria, and on a refined treatment of sub-daily fluctuations of

595

snow depth data. In addition, a mass-conserving procedure has been intro-

596

duced to separate solid and liquid events, and to remedy rain gauges catch

597

deficiency. 26

598

The processing routine and the model need to be coupled, because it is not

599

possible to distinguish between physical variations of SD and noise, at the

600

sub-daily resolution. A processing-modeling routine has been defined, able

601

to filter data and simulate snowpack dynamics at the same time. Multi-year

602

data from 40 sites of the SNOTEL network have been chosen to illustrate

603

the performance of the routine.

604

This application has shown that 1) the processing routine is able to re-

605

move noise, preserving the significant signal in data-series (Section 4.1); 2)

606

the processing-modeling routine succeeds both in correctly predicting snow-

607

pack mass dynamics, and in filtering raw data (Section 4.2); 3) estimated

608

data of nivometric coefficients and catch deficiency are consistent with the

609

literature (Section 4.3.3). An evaluation of parameters’ variability is given

610

in Section 4.3.

611

Processing-modeling low performances (i.e., sites where average N SE

612

values are lower than 0.6) are mainly obtained in the southern and western

613

zones of the study area (Alaska excluded, where only two sites have been

614

considered). In the same locations, ∆TM AX values show a greater variability.

615

The limited amount of stations considered here does not allow inferring any

616

regional pattern of performances of the proposed routine. Future develop-

617

ments should investigate whether the proposed routine performs worse or

618

better based on a regional scale, since this would be of great usefulness to

619

predict a priori the expected performances of the routine.

27

620

6. Acknowledgements

621

Graditude is due to all the SNOTEL technicians involved for the help in

622

R for the data handling. We thank the National Atlas of the United States

623

map used in Figures 1(a) and 1(b). We thank Prof. Jeff Dozier for his

624

helpful comments to the manuscript.

28

625

References

626

[1] Anderson, E.A. 1976. A point Energy and Mass Balance Model of a Snow

627

Cover, NOAA Technical Report NWS 19 , 172 pp.

628

[2] Anderson, J., Wirt, J. 2008. Ultrasonic snow depth sensor accuracy, reli-

629

ability, and performance, Proceedings of the 76th Annual Western Snow

630

Conference, Hood River, Oregon, 99-102.

631

[3] Avanzi, F. 2011. A dynamic model of snowpack density, depth and mass

632

content and its validation with SNOTEL hourly data, M. Sc. Thesis,

633

Politecnico di Milano.

634

[4] Avanzi, F., De Michele, C., Ghezzi, A. 2014. Liquid-Solid Partitioning of

635

Precipitation along an Altitude Gradient and Its Statistical Properties:

636

An Italian Case Study. American Journal of Climate Change, 3, 71-82.

637

DOI:10.4236/ajcc.2014.31007.

638

[5] Avanzi, F., Caruso, M., Jommi, C., De Michele, C., Ghezzi, A. 2014.

639

Continuous-time monitoring of liquid water content in snowpacks using

640

capacitance probes: a preliminary feasibility study. Advances in Water

641

Resources, 3, in press. DOI:10.1016/j.advwatres.2014.02.012.

642

[6] Bales, R.C., Molotch, N.P., Painter, T.H., Dettinger, M.D., Rice, R.,

643

Dozier, J. 2006. Mountain hydrology of the western United States, Water

644

Resources Research, 42(8), W08432. DOI:10.1029/2005WR004387

645

[7] Bardsley, T., Julander, R. 2005. The documentation of extreme events:

646

two case studies in Utah, water year 2005, Proceedings of the 73rd Annual

647

Western Snow Conference, Great Falls, Montana, 33-42. 29

648

[8] Bartelt, P., Lehning, M. 2002. A physical SNOWPACK model for the

649

Swiss avalanche warning Part I: numerical model, Cold Regions Science

650

and Technology, 35, 123-145.

651

[9] Bavay, M., Lehning, M., Jonas, T, L¨owe, H. 2009. Simulations of future

652

snow cover and discharge in Alpine headwater catchments, Hydrological

653

Processes, 23, 95-108. DOI: 10.1002/hyp.7195

654

[10] Bavera, D., Bavay, M., Jonas, T., Lehning, M., De Michele, C. 2014. A

655

comparison between two statistical and a physically-based model in snow

656

water equivalent mapping, Advances in Water Resources, 63, 167178.

657

DOI:http://dx.doi.org/10.1016/j.advwatres.2013.11.011

658

[11] Berghuijs, W. R., Woods, R. A., Hrachowitz, M. 2014. A precipitation

659

shift from snow towards rain leads to a decrease in streamflow, Nature

660

Climate Change. DOI: 10.1038/nclimate2246.

661

[12] Brun, E., David, P., Sudul, M., Brunot, G. 1992. A numerical model to

662

simulate snow-cover stratigraphy for operational avalanche forecasting,

663

Journal of Glaciology, 38(128), 13-22.

664

[13] Carmagnola, C. M., Morin, S., Lafaysse, M., Domine, F., Lesaffre, B.,

665

Lejeune, Y., Picard, G., Arnaud, G. 2014. Implementation and evalu-

666

ation of prognostic representations of the optical diameter of snow in

667

the SURFEX/ISBA-Crocus detailed snowpack model, The Cryosphere,

668

8, 417-437. DOI: 10.5194/tc-8-417-2014

669

[14] Clow, D.W. 2010. Changes in the Timing of Snowmelt and Streamflow

30

670

in Colorado: A Response to Recent Warming, Journal of Climate, 23(9),

671

2293-2306. DOI:10.1175/2009JCLI2951.1

672

[15] Cox, L.M., Bartee, L.D., Crook, A.G., Farnes, P.E., Smith, J.L. 1978.

673

The care and feeding of snow pillows, Proceedings of the 46th Annual

674

Western Snow Conference, Otter Rock, Oregon, 40-47.

675

[16] De Michele, C., Avanzi, F., Ghezzi, A., Jommi, C. 2013. Investigating

676

the dynamics of bulk snow density in dry and wet conditions using a

677

one-dimensional model, The Cryosphere, 7, 433-444. DOI: 10.5194/tc-7-

678

433-2013

679

[17] De Michele, C., Ignaccolo, M. 2013. New perspectives on rain-

680

fall from a discrete view, Hydrological Processes, 27(16), 23792382.

681

DOI:10.1002/hyp.9782.

682

683

[18] Duchon, C., Biddle, C. 2010. Undercatch of tipping-bucket gauges in high rain rate events, Advances in Geosciences, 25, 11-15.

684

[19] Egli, L., Jonas, T., Meister, R. 2009. Comparison of different automatic

685

methods for estimating snow water equivalent, Cold Regions Science and

686

Technology, 57, 107-115. DOI:10.1016/j.coldregions.2009.02.008

687

[20] Fassnacht, S.R., Dressler, K.A., Bales, R.C. 2003. Snow water

688

equivalent interpolation for the Colorado River Basin from snow

689

telemetry (SNOTEL) data, Water Resources Research, 39(8), 1208.

690

DOI:10.1029/2002WR001512.

691

[21] Fassnacht, S.R. 2004. Estimating Alter-shielded gauge snowfall un-

692

dercatch, snowpack sublimation, and blowing snow transport at six 31

693

sites in the coterminous USA, Hydrological Processes, 18, 34813492.

694

DOI:10.1002/hyp.5806.

695

[22] Fassnacht, S.R. 2006. Upper versus lower Colorado River sub-basin

696

streamflow: characteristics, runoff estimation and model simulation, Hy-

697

drological Processes, 20(10), 2187-2205. DOI:10.1002/hyp.6202.

698

[23] Fassnacht, S.R., Derry, J.E. 2010. Defining Similar Regions of Snow in

699

the Colorado River Basin using Self Organizing Maps (SOMs). Water

700

Resources Research, 46, W04507. DOI:10.1029/2008WR007835.

701

[24] Groisman, P. Y., Legates, D. R. 1994. The accuracy of United States pre-

702

cipitation data, Bulletin of the American Meteorological Society, 75(3),

703

215-227.

704

[25] Gubler, H. 1981. An inexpensive remote snow-depth gauge based on

705

ultrasonic wave reflection from the snow surface, Journal of Glaciology,

706

27(95), 157 163.

707

[26] Haberlandt, U. 2007. Geostatistical interpolation of hourly precipita-

708

tion from rain gauges and radar for a large-scale extreme rainfall event,

709

Journal of Hydrology, 332, 144157. DOI:10.1016/j.jhydrol.2006.06.028.

710

[27] Hamlet, A.F., Mote, P.W., Clark, M.P., Lettenmaier, D.P. 2005. Ef-

711

fects of Temperature and Precipitation Variability on Snowpack Trends

712

in the Western United States, Journal of Climate, 18(21), 45454561.

713

DOI:10.1175/JCLI3538.1 .

714

[28] Hock, R. 2003. Temperature index melt modelling in mountain areas,

715

Journal of Hydrology, 282, 104115. DOI:10.1016/S0022-1694(03)00257-9. 32

716

[29] Johnson, J.B., Schaefer, G.L. 2002. The influence of thermal, hy-

717

drologic, and snow deformation mechanisms on snow water equiva-

718

lent pressure sensor accuracy, Hydrological Processes, 16(18), 35293542.

719

DOI:10.1002/hyp.1236

720

721

[30] Johnson, J.B. 2004. A theory of pressure sensor performance in snow, Hydrological Processes, 18(1), 5364. DOI:10.1002/hyp.1310

722

[31] Johnson, J.B., Marks D. 2004. The detection and correction of snow

723

water equivalent pressure sensor errors, Hydrological Processes, 18(18),

724

35133525. DOI:10.1002/hyp.5795

725

[32] Jonas, T., Marty, C., Magnusson, J. 2009. Estimating the snow water

726

equivalent from snow depth measurements in the Swiss Alps, Journal of

727

Hydrology, 378, 161-167. DOI:10.1016/j.jhydrol.2009.09.021

728

[33] Jordan, R. 1991. A one-dimensional temperature model for a snow cover,

729

Technical documentation for SNTHERM.89, Spec. Rep. 91-16, US Army

730

Corps of Eng., Cold Regions Res. Eng. Lab., Hanover.

731

[34] Julander, R.P. 2007. Soil Surface Temperature Difference Between Steel

732

and Hypalon Pillows, Proceedings of the 75th Annual Western Snow Con-

733

ference, Kailua-Kona, Hawaii.

734

[35] Julander, R.P., Wilson, G.R., Nault, R. 1998. The Franklin basin prob-

735

lem, Proceedings of the 66th Annual Western Snow Conference, Snow-

736

bird, Utah, 153-156.

737

[36] Katsushima, T., Kumakura, T., Takeuchi, Y. 2009. A multiple snow 33

738

layer model including a parameterization of vertical water channel pro-

739

cess in snowpack, Cold Regions Science and Technology, 59, 143-151.

740

DOI:10.1016/j.coldregions.2009.09.002

741

[37] Kienzle, S.W. 2008. A new temperature based method to separate rain

742

and snow, Hydrological Processes, 25, 5067-5085. DOI:10.1002/hyp.7131

743

[38] Kondo, J., Yamazaki, T. 1990. A prediction model for snowmelt, snow

744

surface temperature and freezing depth using a heat balance method,

745

Journal of Applied Meteorology, 29, 375-384.

746

[39] Kumar, M., Marks, D., Dozier, J., Reba, M., Winstral, A. 2013.

747

Evaluation of distributed hydrologic impacts of temperature-index and

748

energy-based snow models, Advances in Water Resources, 56, 77-89.

749

DOI:10.1016/j.advwatres.2013.03.006

750

[40] Larson, L.W., Peck, E.L. 1974. Accuracy of Precipitation Measurements

751

for Hydrologic Modeling, Water Resources Research, 10(4), 857-863.

752

[41] Lehning, M., Bartelt, P., Brown, B., Fierz, C., Satyawali, P. 2002. A

753

physical SNOWPACK model for the Swiss avalanche warning Part II:

754

Snow microstructure, Cold Regions Science and Technology, 35, 147-167.

755

[42] Lehning, M., Bartelt, P., Brown, B., Fierz, C. 2002. A physical SNOW-

756

PACK model for the Swiss avalanche warning Part III: meteorological

757

forcing, thin layer formation and evaluation, Cold Regions Science and

758

Technology, 35, 169-184.

759

[43] Leibowitz, S. G., Wigington Jr., P. J., Comeleo, R. L., Ebersole, J. L. 34

760

2012. A temperature-precipitation-based model of thirty-year mean snow-

761

pack accumulation and melt in Oregon, USA, Hydrological Processes,

762

26(5), 741759. DOI: 10.1002/hyp.8176

763

[44] Lenderink, G., van Meijgaard, E. 2008. Increase in hourly precipitation

764

extremes beyond expectations from temperature changes, Nature Geo-

765

science, 1(8), 511-514. DOI:10.1038/ngeo262

766

[45] Luce, C.H., Tarboton, D.G., Cooley, K.R. 1998. The influence of the

767

spatial distribution of snow on basin-averaged snowmelt, Hydrological

768

Processes, 12, 1671-1683.

769

[46] Luce, C.H., Tarboton, D.G., Cooley, K.R. 1999. Sub-grid parameteri-

770

zation of snow distribution for an energy and mass balance snow cover

771

model, Hydrological Processes, 13, 1921-1933.

772

[47] Macelloni, G., Paloscia, S., Pampaloni, P., Brogioni, M., Ranzi, R.,

773

Crepaz, A. 2005. Monitoring of Melting Refreezing Cycles of Snow With

774

Microwave Radiometers: The Microwave Alpine Snow Melting Exper-

775

iment (MASMEx 2002 - 2003), IEEE Transactions on Geoscience and

776

Remote Sensing, 43(11), 2431-2442.

777

[48] McCuen, R. H., Knight, Z., Cutter, A. G. 2006. Evaluation of the Nash-

778

Sutcliffe Efficiency Index, Journal of Hydrologic Engineering, 11, 597-

779

602. DOI:10.1061/(ASCE)1084-0699(2006)11:6(597)

780

[49] McGinnis, D.L. 1997. Estimating Climate-Change Impacts on Colorado

781

Plateau Snowpack Using Downscaling Methods, Professional Geographer ,

782

49, 117-125. 35

783

[50] McGurk, B.J. 1986. Precipitation and snow water equivalent sensors: an

784

evaluation, Proceedings of the 54th Annual Western Snow Conference,

785

Phoenix, Arizona, 71-80.

786

[51] Meromy, L., Molotch, N.P., Link, T.E., Fassnacht, S.R., Rice, R. 2012.

787

Subgrid variability of snow water equivalent at operational snow stations

788

in the western USA, Hydrological Processes. DOI:10.1002/hyp.9355

789

[52] Mizukami,

N.,

Perica,

S. 2008. Spatiotemporal Characteristics

790

of Snowpack Density in the Mountainous Regions of the West-

791

ern United States, Journal of Hydrometeorology, 9(6), 1416-1426.

792

DOI:10.1175/2008JHM981.1

793

[53] Molotch, N. P., Fassnacht, S. R., Bales, R. C., Helfrich, S. R. 2004.

794

Estimating the distribution of snow water equivalent and snow extent

795

beneath cloud cover in the SaltVerde River basin, Arizona, Hydrological

796

Processes, 18(9), 1595-1611. DOI:10.1002/hyp.1408

797

[54] Morin, S., Lejeune, Y., Lesaffre, B., Panel, J.-M., Poncet, D., David,

798

P., Sudul, M. 2012. An 18-yr long (1993-2011) snow and meteorological

799

dataset from a mid-altitude mountain site (Col de Porte, France, 1325m

800

alt.) for driving and evaluating snowpack models, Earth Syst. Sci. Data,

801

4, 13-21. DOI:10.5194/essd-4-13-2012

802

[55] Mote, P.W. 2003. Trends in snow water equivalent in the Pacific North-

803

west and their climatic causes, Geophysical Research Letters, 30(12),

804

1601. DOI:10.1029/2003GL017258

36

805

[56] Nash, J.E.., Sutcliffe, J.V. 1970. River flow forecasting through concep-

806

tual models. Part I: a discussion of principles, Journal of Hydrology, 10,

807

282-290.

808

[57] Nomura, M. 1994. Studies on the Delay Mechanism of Runoff to

809

Snowmelt, Contributions from the Institute of Low Temperature Science,

810

Hokkaido University, 39, 1-49.

811

[58] Osterhuber, R.S., Edens, T., McGurk, B.J. 1994. Snow depth measure-

812

ment using ultrasonic sensors and temperature correction, Proceedings

813

of the 62nd Annual Western Snow Conference, Santa Fe, New Mexico,

814

159-162.

815

[59] Pfeffer, W.T., Illangasekare, T.H., Meier, M.F. 1990. Analysis and

816

modeling of melt-water refreezing in dry snow, Journal of Glaciology,

817

36(123).

818

[60] Raleigh, M.S., Lundquist, J.L. 2012. Comparing and combining SWE

819

estimates from the SNOW-17 model using PRISM and SWE reconstruc-

820

tion, Water Resources Research, 48(1). DOI:10.1029/2011WR010542

821

[61] Richards, R.P. 1983. Snow pillow overloading by slope induced creep,

822

Proceedings of the 52nd Annual Western Snow Conference, Sun Valley,

823

Idaho, 194-187.

824

[62] Ryan, W.A., Doesken, N.J., Fassnacht, S.R. 2008. Preliminary results of

825

ultrasonic snow depth sensor testing for National Weather Service (NWS)

826

snow measurements in the US, Journal of atmospheric and oceanic tech-

827

nology, 25, 2748-2757. 37

828

[63] Sade, R., Rimmer, A., Iggy Litaor, M., Shamir, E., Furman A.

829

2011. The sensitivity of snow-surface temperature equation to sloped

830

terrain, Technical Note, Journal of Hydrology, 408, 308-313. DOI:

831

10.1016/j.jhydrol.2011.08.001

832

833

834

[64] Schaefli, B., Gupta, H. V. 2007. Do Nash values have value?, Hydrological Processes, 21, 2075-2080. DOI: 10.1002/hyp.6825 [65] Schleef,

S.,

L¨owe,

H.

2013.

X-ray

microtomography

analy-

835

sis of isothermal densification of new snow under external me-

836

chanical stress,

837

http://dx.doi.org/10.3189/2013JoG12J076.

Journal of Glaciology,

59(214),

233-243. DOI:

838

[66] Schmid, L., Heilig, A., Mitterer, C., Schweizer, J., Maurer, H., Okorn,

839

R., Eisen, O. 2014. Continuous snowpack monitoring using upward-

840

looking ground-penetrating radar technology, Journal of Glaciology, 60,

841

509-525. DOI: 10.3189/2014JoG13J084.

842

[67] Schmid, M.O., Gubler, S., Fiddes, J., Gruber, S. 2012. Inferring

843

snowpack ripening and melt-out from distributed measurements of

844

near-surface ground temperatures, The Cryosphere, 6(5), 1127-1139.

845

DOI:10.5194/tc-6-1127-2012

846

[68] Schneebeli, M. 1995. Development and stability of preferential flow paths

847

in a layered snowpack, Biogeochemistry of seasonally snow-covered catch-

848

ments (Proceedings of a Boulder Symposium, July 1995), IAHS publ. no.

849

228.

38

850

[69] Serreze, M.C., Clark, M.P., Armstrong, R.L., McGinnis, D.A., Pulwarty,

851

R.S. 1999. Characteristics of the western United States snowpack from

852

snowpack telemetry (SNOTEL) data, Water Resources Research, 35(7),

853

2145-2160. DOI:10.1029/1999WR900090

854

[70] Smith, F.W., Boyne, H.S. 1981. Snow pillow behavior under controlled

855

laboratory conditions, Proceedings of the 49th Annual Western Snow

856

Conference, St. George, Utah, 13-22.

857

[71] Tanner, B.D., Gaza, B. 1990. Automated snow depth and snowpack tem-

858

perature measurements, Proceedings of the 58th Annual Western Snow

859

Conference, Sacramento, California. 73-78.

860

[72] Tarboton, D. G., Luce, C. H. 1996. Utah energy balance snow accumu-

861

lation and melt model (UEB), Computer model technical description and

862

users guide, Utah Water Res. Lab., Logan.

863

[73] Tobin, C., Schaefli, B., Nicotina, L., Simoni, S., Barrenetxea,

864

G., Smith, R., Parlange, M., Rinaldo, A. 2013. Improving the

865

degree-day method for sub-daily melt simulations with physically-

866

based diurnal variations, Advances in Water Resources, 55, 149-164.

867

DOI:http://dx.doi.org/10.1016/j.advwatres.2012.08.008.

868

[74] Vionnet, V., Brun, E., Morin, S., Boone, A., Faroux, S., Le Moigne, P.,

869

Martin, E., Willemet, J.-M. 2012. The detailed snowpack scheme Crocus

870

and its implementation in SURFEX v7.2, Geosci. Model Dev., 5, 773791.

871

DOI:10.5194/gmd-5-773-2012.

39

872

[75] Wever, N., Fierz, C., Mitterer, C., Hirashima, H., Lehning, M. 2014.

873

Solving Richards Equation for snow improves snowpack meltwater runoff

874

estimations in detailed multi-layer snowpack model, The Chryosphere, 8,

875

257-274. DOI:10.5194/tc-8-257-2014.

876

[76] Yamaguchi, S., Nakai, S., Iwamoto, K., Sato, A. 2009. Influence of

877

Anomalous Warmer Winter on Statistics of Measured Winter Precipi-

878

tation Data, Journal of Applied Meteorology and Climatology, 48, 2403-

879

2409.

880

[77] Zhang, Y., Wang, S., Barr, A.G., Black, T.A. 2008. Impact of

881

snow cover on soil temperature and its simulation in a boreal as-

882

pen fores, Cold Regions Science and Technology, 52(3), 355-370.

883

DOI:10.1016/j.coldregions.2007.07.001

40

Table 1: The selected SNOTEL sites. ID

SNOTEL number

SNOTEL name

State

Elevation

Dataset

S1

301

Adin Mtn

CAL

1886

2005-2011

S2

587

Lobdell Lake

CAL

2814

2005-2011

S3

539

Ind. Camp

CAL

2134

2002-2011

S4

724

Rubicon #2

CAL

2343

2004-2011

S5

705

Promontory

ARI

2417

2005-2011

S6

308

Baker Butte

ARI

2200

2003-2011

S7

854

Wesner Springs

NEWM

3389

2005-2009

S8

394

Chamita

NEWM

2560

2006-2011

S9

922

Santa Fe

NEWM

3488

2008-2011

S10

498

Granite Peak

NEV

2600

2005-2011

S11

476

Fawn Creek

NEV

2133

2008-2011

S12

849

Ward Mountain

NEV

2800

2008-2011

S13

766

Snowbird

UTAH

2938

2004-2011

S14

631

Mining Fork

UTAH

2500

1999-2011

S15

844

Vernon Creek

UTAH

2255

2001 - 2011

S16

983

Clayton Springs

UTAH

3062

2001-2011

S17

378

Burro Mountain

COL

2865

2004 - 2011

S18

580

Lily Pond

COL

3350

2004-2011

S19

467

Elk River

COL

2650

2004-2011

S20

341

Big Red Mountain ORE

1844

2005-2011

S21

800

Summer Rim

2157

2004-2011

ORE

Table 1: The table continues in the next page

41

Table 1: The table continues from the previous page

ID

SNOTEL number

SNOTEL name

State

Elevation

Dataset

S22

653

Mt. Howard

ORE

2410

2004-2011

S23

655

Mud Ridge

ORE

1240

2007-2010

S24

644

Moses Mtn

WAS

1527

2004-2011

S25

502

Green Lake

WAS

1804

2005-2011

S26

817

Thunder Basin

WAS

1316

2008-2011

S27

988

Hidden Lake

IDA

1536

2003-2011

S28

411

Cool Creek

IDA

1914

2006-2010

S29

490

Galena Summit

IDA

2676

2007-2011

S30

830

Trinity Mountain

IDA

2368

2006-2009

S31

903

Nevada Ridge

MON

2139

2006-2011

S32

783

Sleeping Woman

MON

1874

2006-2011

S33

347

Black Bear

MON

2490

2006-2011

S34

872

Windy Peak

WYO

2407

2008-2011

S35

822

Togwotee Pass

WYO

2919

2004-2010

S36

499

Grassy Lake

WYO

2214

2007-2011

S37

367

Brooklin Lake

WYO

3100

2007-2011

S38

354

Blind Park

SDAK

2100

2007-2011

S39

966

Kenai Moose Pens

ALA

91

2007-2011

S40

1096

May Creek

ALA

490

2008-2010

884

42

Table 2: Model variables and parameters.

Symbol

Definition

hS

Depth of the ice structure (m)

ρD

Snowpack dry density (kg/m3 )

hW

Depth of the liquid constituent of the snowpack (m)

ρN S

Density of new snow (kg/m3 , [1])

s

New solid events (m/h), input data

I (T, hS )

Function equal to 1 if TA > 0 ◦ C, which tends to zero if hS → 0

T

Air temperature (◦ C), input data

Tτ

Threshold temperature of melting (0◦ C)

p

New liquid events (m/h), input data

ρW

Liquid water density (1000 kg/m3 )

θ

Mean volumetric liquid water content

d

1.25, empirical coefficient [57]

c1

0.01 m2 /h/kg, empirical coefficient [77]

TS

Mean snow temperature, retrieved from air temperature [38]

n

Porosity of ice matrix

hi

Maculay’s brackets

h = hS + hhW − nhS i

Control depth (m)

ρ=

ρD hS +ρW hW h

SW E =

ρh ρW

Bulk snow density (kg/m3 ) Snow water equivalent (mm)

a

Degree hour parameter (m/h/◦ C)

c

Runoff parameter (1/md−1 )/h

43

885

Table 3: Model equations. The symbols reported in the equations are given in Table 2.

Equations dhS dt

= − ρhDS ρdtD + ρD ρW

ρN S s ρD

− I (T, hS ) (a (T − Tτ ))

dhW dt

=p+

dρD dt

= c1 ρ2D hS e(0.08(T −Tτ )−0.021ρD ) +

I (T, hS ) (a (T − Tτ )) − cθhdW

886

44

(ρN S −ρD ) s hS

(a)

(b) Figure 1: A map of western United States, panel (a), and of Alaska, panel (b), with the location of the SNOTEL sites used in this analysis.

45

Figure 2: Elevation distribution of the 40 SNOTEL sites selected in this analysis. The two sites, at an elevation lower than 800 m a.s.l., are situated in Alaska.

46

Figure 3: The proposed processing routine for the SNOTEL hourly data of snow depth, snow water equivalent, air temperature and cumulative precipitation.

47





 
































 



 

 



 










(b)



















(a)





 










          

 







 



(c)

(d)








 





 



(e) Figure 4: An example of SNOTEL hourly data, for site S26 and W. Y. 2011: Panel a) shows raw (green) versus processed (black) snow depth , Panel b) gives raw (green) versus processed (black) cumulative precipitation, solid (purple) and liquid (blue) cumulative precipitation , Panel c) reports raw (green) versus processed (black) data of SW E, Panel d) provides a zoom of panel a), while Panel e) reports raw (green) versus processed (black) data of hourly positive differences in cumulative precipitation.

48 .









 




























(a)
 








 













(b)








 








 

(c) Figure 5: An example of multi-year dynamics of snow depth, bulk snow density and SWE, at site S33. Panel a) gives processed SD data (black), modeled h (red) and modeled hW (blue), Panel b) shows processed bulk snow density data (black), modeled bulk snow density (red), and modeled dry snow density (blue), Panel c) provides processed (black) and modeled (red) SW E.

49






















 










 


















(b)



 






 







(a)






 














 













(c)

(d)













 









 








 





 

(e)

(f)

Figure 6: Comparison between modeled results obtained using directly raw data as input (panel a, c and e) versus processing-modeling results (panel b, d and f) at site S26 and W.Y. 2009. Colors are the same as in Figure 5. In Panel (a), (c) and (e), processed data are substituted with raw data (green). In the same panels, they are reported in purple modeled results one would obtain processing raw data-series using visual thresholds only, i.e. not considering ∆TM AX and MLAST .

50













  

 




 





 



Figure 7: Box-plot of the Nash-Sutcliffe coefficients NSE of the processing-modeling routine using the first-year calibration procedure, over all the 40 sites. On the left side, NSESD , NSEρ and NSESWE , for each site, have been calculated at the calibration year. On the right side, NSESD , NSEρ and NSESWE are mean values calculated on the whole dataset of each site.

51

(a)

(b)

(c) Figure 8: Box-plots of the average coefficient of determination NSE of the processingmodeling routine using the one-year calibration procedure. NSESD in panel a), NSEρ in panel b), and NSESWE in panel c).

52

Figure 9: An example of ∆TM AX and a calibration. Panel a) gives the error surface, calculated at site S4, for the water year 2004. Panel b) shows contours of the same error surface. a values have been multiplied by -1 for sake of clarity.

53

(a)

(b)

(c) Figure 10: Box-plot of the a parameter, panel (a), of the ∆TM AX parameter, panel (b), and the c parameter, panel (c), for the 40 sites.

54







  




































   (a)
























 
  (b)

Figure 11: Panel (a): Median values of the a parameter plotted as a function of the site elevation. Panel (b): Values of the c parameter plotted as a function of the site elevation.

55






















 
 (a)
















          
 (b) Figure 12: Nivometric coefficient of each calibration year plotted as a function of the site elevation, panel (a), and ratio between PM OD and observed P precipitation, of each calibration year, plotted as a function of the nivometric coefficient of the same calibration year.

56

,ŝŐŚůŝŐŚƚƐ͗ • • • • •

tĞĞǀĂůƵĂƚĞƚŚĞƵƐĞŽĨ^EKd>ŚŽƵƌůLJĚĂƚĂŝŶĂƐŶŽǁĚLJŶĂŵŝĐŵŽĚĞů͖ ƵĞƚŽĚĂƚĂŶŽŝƐĞ͕ĂƉƌŽĐĞƐƐŝŶŐƌŽƵƚŝŶĞŚĂƐďĞĞŶĚĞƐŝŐŶĞĚ͖ dŚĞƌŽƵƚŝŶĞŝƐĐŽƵƉůĞĚƚŽƚŚĞŵŽĚĞůƚŽƌĞŵŽǀĞĞƌƌŽŶĞŽƵƐŚŽƵƌůLJŽƐĐŝůůĂƚŝŽŶƐ͖ ϰϬ^EKd>ƐŝƚĞƐƚŚƌŽƵŐŚŽƵƚǁĞƐƚĞƌŶh͘^͘ŚĂǀĞďĞĞŶƐĞůĞĐƚĞĚĂƐĂƚĞƐƚ͖ WĞƌĨŽƌŵĂŶĐĞƐŽĨƚŚĞƌŽƵƚŝŶĞĂƌĞŐŽŽĚ͕ďŽƚŚŝŶƉƌŽĐĞƐƐŝŶŐĂŶĚŝŶŵŽĚĞůŝŶŐ͘