Bulk volumetric liquid water content in a seasonal snowpack: modeling its dynamics in different climatic conditions

Advances in Water Resources 86 (2015) 1–13

Contents lists available at ScienceDirect

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

Bulk volumetric liquid water content in a seasonal snowpack: modeling its dynamics in different climatic conditions Francesco Avanzi a,∗, Satoru Yamaguchi b, Hiroyuki Hirashima b, Carlo De Michele a a b

Department of Civil and Environmental Engineering, Politecnico di Milano, Piazza Leonardo da Vinci, Milano, Italy Snow and Ice Research Center, National Research Institute for Earth Science and Disaster Prevention, Suyoshi-machi, Nagaoka-shi, Niigata-ken, 940-0821, Japan

a r t i c l e

i n f o

Article history: Received 3 April 2015 Revised 17 September 2015 Accepted 23 September 2015 Available online 5 October 2015 Keywords: Liquid water content Degree-hour Temperature-index Refreezing Snow Numerical model

a b s t r a c t We focus on the dynamics of volumetric liquid water content in seasonal snow covers. This is a key variable describing the fate of snowpacks during the melting season. However, its measurement and/or prediction by means of models at high spatial and temporal resolutions is still difficult due to both practical and theoretical reasons. To overcome these limitations in operational applications, we test the capability of a onedimensional model to predict the dynamics of bulk volumetric liquid water content during a snow season. Multi-year data collected in three experimental sites in Japan are used as an evaluation. These sites are subjected to different climatic conditions. The model requires the calibration of one or two parameters, according to the degree of detail used. Either a simple temperature-index or a coupled melt–freeze temperature-index approach are considered to predict melting and/or melt–freeze dynamics of liquid water. Results show that, if melt–freeze dynamics are modeled, median absolute differences between data and predictions are consistently lower than 1 vol% at the sites where data of liquid water content are available. In addition, we find also that the model predicts correctly a dry condition in 80% of the observed cases at a site where calibration data are scarce. At the same site, observed isothermal conditions of the snow cover at 0 °C correspond to predictions of bulk volumetric liquid water content that are greater than 0. © 2015 Elsevier Ltd. All rights reserved.

1. Introduction Snow is a multiphase mixture of three main constituents: ice, air and liquid water [16]. Ice, in the form of snow grains, constitutes a solid skeleton with high porosities [44,66], while liquid water occupies a certain amount of pores volume. Following mixture theory [58], volumetric liquid water content (henceforth LWC, or θ W ) is nowadays widely used as a measure of liquid water amount in a snow cover [2]. This is defined as the ratio between the volume occupied by liquid water and the total volume of the mixture, at a given point (vol%) [24]. Liquid water dynamics in snow have major impacts on wet snow avalanche prediction [4,54,56,73], runoff modeling [80] or snow albedo calculation [75]. All these applications benefit from precise estimations of LWC. However, continuous-time, automatic and undisturbed measurements of LWC are still difficult [2,73]. Consequently, rare examples exist of this type of measurement (see as an example the works by [28,29,45,55,67]).

∗

Corresponding author. Tel.: +39223996235. E-mail address: [email protected] (F. Avanzi).

http://dx.doi.org/10.1016/j.advwatres.2015.09.021 0309-1708/© 2015 Elsevier Ltd. All rights reserved.

On the other hand, the modeling of LWC dynamics in snow is notoriously challenging [80]. This is motivated by the need to couple hydraulic [8,70,86], mechanical [52,74] and thermal [35,62,76] processes. Moreover, it is well-known that some components of liquid water routing in snow develop within a multi-dimensional domain [32]. In fact, liquid water percolation in seasonal snow is characterized by the development of preferential flow paths [39,40,51,68,78], lateral or slope-parallel flow [19], the development of capillary barriers [12,31,56,78,80,83,84], and the occurrence of melt–freeze dynamics [62]. Moreover, snow structure changes with time [18,37,63,69], and liquid water induces additional quick modifications of grain shape and size. This peculiar process calls wet snow metamorphism [6,11]. As a result, a comprehensive modeling of LWC dynamics along a vertical snow profile is still an open issue for many snow models due to both practical and theoretical reasons. We consider here a simple model that explicitly estimates θ W , namely the bulk volumetric liquid water content of the snow cover at a given point. Input and calibration data requirements are relatively low and this may provide a practical tool to obtain predictions over wide areas where input data are restricted. Considering volumetric LWC in addition to liquid water mass (which is already included by, e.g., some distributed snow hydrology models for runoff prediction [25,65]) provides a quantity that can be directly used in

2

F. Avanzi et al. / Advances in Water Resources 86 (2015) 1–13 Table 2 Details about the instrumentation used at each site (data are all at an hourly resolution). T is air temperature, P is total precipitation, W is wind speed, SD is snow depth, SWE is Snow Water Equivalent, R is snowmelt runoff, while LWC is volumetric liquid water content. Variable

Nagaoka

Shinjo

Sapporo

T P

Platinum therm. Tipping bucket rain gauge Aerovane Light wave phase difference method Pressure pillow Lysimeter Denoth meter

Platinum therm. Tipping bucket rain gauge Aerovane Light wave phase difference method Pressure pillow Lysimeter Denoth meter

Thermo–hygrom. Tipping bucket rain gauge Aerovane Ultrasonic depth sensor – – –

W SD SWE R LWC

Fig. 1. Location of the three study sites considered. Table 1 Sites location, type of snow region according to the classification by [36], water years (WY) considered in model calibration, and water years considered in θ W evaluation. Sites name

Coordinates

Snow region

WY (cal.)

WY (θ W eval.)

Nagaoka

37° 25’ N 138° 53’ E 97 m ASL 38° 47’ N 140° 19’ E 127 m ASL 43° 5’ N 141° 20’ E 15 m ASL

Wet

1997–2004 2006–2014

Intermediate

2006 2008–2013

1997/2001–2004 2006 2008–2014 2006 2008–2013

Dry

2006–2013

2006–2013

Shinjo

Sapporo

wet snow avalanches forecasting [57]. Moreover, θ W is frequently used to classify the wetness condition of a snow cover [24,73] and to describe the transition between flow regimes [14]. In the model, either a simple temperature-index approach or a combined melt– freeze temperature-index approach are included to predict the melting and/or melt–freeze mass fluxes. As an evaluation, we compare these estimations with measurements of the vertical properties of a seasonal snow cover in three different sites in Japan, where multiyear data to calibrate, force and evaluate the proposed model are available. These are subjected to different climatic conditions with respect to θ W dynamics. 2. Data To evaluate model performances, we consider here three experimental sites in Japan. These are the Snow and Ice Research Center in Nagaoka (henceforth Nagaoka [81]), the Shinjo Cryospheric Environment Laboratory, Snow and Ice Research Center, NIED in Shinjo (henceforth Shinjo, see e.g. [59]) and the Institute of Low Temperature Science Hokkaido University in Sapporo (henceforth Sapporo [1,30]). These stations are all placed in flat grasslands that are characterized by a regular topography. Recently, [33] used data from these sites to evaluate SNOWPACK performances in Japan. Fig. 1 reports the location of the three sites, while Table 1 reports details about the water years considered in the analysis. We define a water year as a period that starts on October 1 and ends on following September 30. Here, we will indicate a water year using the year when it ends (e.g., the water year 2006/2007 is indicated as 2007). The same table reports the classification of each site according to the climatic division of snowy areas in Japan by [36]. Accordingly, Nagaoka is located in the wet-snow region. This means that “every layer

of deposited snow is wet throughout the winter season” (see Table 2 in [36]), and that monthly average temperature during January (i.e., the coldest month during the year) is positive. This represents an unusual climatic condition with respect to alpine areas [71]. On the contrary, Sapporo falls in the dry-snow region, where “deposited snow is dry at least in mid-winter” (see again Table 2 in [36]), while Shinjo belongs to an intermediate snow region between them. These three sites represent an interesting case study to evaluate the performances of a simple model in evaluating LWC in different climatic conditions. At each site, hourly data of wind speed W, total precipitation P, air temperature T and snow depth SD are available during each snow season (i.e, usually between December and April, both limits included). In Nagaoka and Shinjo, continuous-time data of Snow Water Equivalent SWE (hence, bulk density ρ ) are also measured during the same periods. Hourly snow runoff R measurements by a snow lysimeter [85] are available during the period 2008–2014 in Nagaoka (9 m2 area), and for the water years 2012 and 2013 in Shinjo (4 m2 area), again during snow seasons. We report in Table 2 an inventory of the instruments used at each site to monitor all the mentioned variables. Regular snow pit measurements are performed at each site during the snow season. Temporal intervals between these measurements are usually between 1 and 10 days. These provide a vertical description of the properties of the snow cover (layers thickness, density, SWE, temperature and grain type, see [24]). In Nagaoka and Shinjo, quantitative measurements of LWC (percentage by mass W (z)) at different elevations z are also collected during snow pit operations. Measurements are usually made using a Denoth meter [15]. This instrument has been widely used to measure LWC profiles in snow [73]. As a rule of thumb, one measurement every 5 or 10 cm is usually collected, but the vertical resolution of these measurements varies depending on layering and simultaneous snow temperature. We have converted this quantity to θ W by evaluating a weighted average of W (z) as a function of measurement distances (W ), and then by using the observed bulk snow density ρ as follows [16]:

θW = W ·

ρ . ρW

(1)

ρ W is the liquid water density (equal to 1000 kg/m3 ). No continuous-time measurement of LWC is available. At each site, pit measurements are performed on the same plot of the meteorological station. In Nagaoka, the distance between the meteorological station and pit locations is approximately 10 m (northern east direction from the meteorological station). In Shinjo, the meteorological station and pit locations are approximately 20 m far (northern east direction from the meteorological station). In Sapporo, this distance is approximately 50 m (southern direction). Therefore, the source of meteorological and pits data are close to each other and characterized by the same top ground surface. Rain-gauge deficiency has been corrected at each location by accounting for wind speed effect. We used the approach suggested in

F. Avanzi et al. / Advances in Water Resources 86 (2015) 1–13

[82] for this purpose. In Nagaoka, a cross-check with SWE increments was used as an additional correction during the water year 2010, due to severe under-catch. Then, events have been classified as solid or liquid using a static threshold approach [43]. Accordingly, we consider an event as rain when simultaneous air temperature is greater than a site-specific threshold (+0.5 °C in Nagaoka and Shinjo, due to past studies [82], and 0 °C in Sapporo), and as snow otherwise. In Sapporo, due to abundant instrumental under-catch, precipitation data were substituted with simultaneous measurements at a close station managed by the Japan Meteorological Agency, and then corrected. This station is located 4 km away from Hokkaido University (southern direction). 3. Methods In this Section, we firstly introduce the model, and, secondly, we discuss its calibration and evaluation in Japan. We discuss both the melt temperature-index approach (HyS model) and the coupled melt–freeze temperature-index approach (HyS 2.0 model). This model is a development of the one originally formulated by [13]. It was evaluated in western US by [3] and recently applied by [26] in Switzerland. In Table 3, all symbols used are given. 3.1. The model Let us consider the snowpack at a given point as an isolated onelayer mixture of a dry and a wet constituent. We define as MS the mass of the dry component (snow grains), and as MW the mass of the wet constituent (liquid water). At any given time, the snow mixture occupies a total volume equal to V, with unitary area and height equal to h. The dry constituent occupies a solid deformable skeleton, composed by snow grains and pores, with total volume equal to VS , unitary area and height equal to hS . Parallely, liquid water occupies a volume equal to VW , unitary area and height equal to hW . We define ρD = MS /VS as the dry density of the snowpack, while ρW = MW /VW = 1000 kg/m3 as the density of water. During most of the season, hW ≤ n · hS , where n = VP /VS = 1 − ρD /ρi is the porosity, VP the total volume of pores in the snowpack, and ρ i is the ice density, equal to 917 kg/m3 . During the last instants of a snow season, due to gradual phase conversion from ice to water, liquid water can be the prevalent phase, and hW > n · hS . We therefore define h = hS + hW − nhS , where  are Macaulay brackets, and provide the argument if this is positive, otherwise 0. Consequently, bulk snow density ρ = (MS + MW )/V = (ρD · hS + ρW · hW )/h, SW E = ρ · h/ρW [16], and θW = hW /h (bulk volumetric liquid water content of the mixture). 3.1.1. Melt temperature-index approach (HyS model) We firstly give the mass balance equations in the integral form (namely, referred to the column of snow) neglecting refreezing. These are as follows (t is time, in hours):

dMS =S−M dt

(2)

dMW =L+M−R dt

(3)

where S is the solid precipitation mass flux, M is the melting mass flux, L the liquid precipitation mass flux, and R the runoff term. All fluxes are in kg/m2 /h. Due to the lack of simplified methods for their computation, these equations do not consider sublimation, evaporation and wind redistribution (see [13] for details). We define S as ρ NS · s, where s is snow precipitation rate (in m/h), while ρNS = 3.6 · W − 0.2 · T + 62 is new snow density, in kg/m3 , following the formula proposed by [38,64] for Japanese areas. W is wind speed in m/s and T is air temperature in °C. L = ρW · p, where

3

Table 3 List of all the symbols used. Symbol

Definition

a

Degree hour parameter (a¯ is its median value) Constants (see Eq. (4)) Modulus of the difference between observed and simulated LWC melt–freeze factor (e¯ is its median value) Snow viscosity Parameters of the viscosity law [77] Melt–freeze mass flux Total height of the snow mixture Height of melt–freeze ice Height of the porous matrix Height of the liquid constituent of the snowpack Multiplicative functions (see Section 3) Constant (see Section 3.1.2) Intrinsic permeability of snow Intrinsic permeability of water in snow Liquid precipitation mass flux Bulk volumetric liquid water content Melting mass flux Mass of melt–freeze ice Mass of the dry constituent of the snowpack Mass of the wet constituent of the snowpack Porosity of the ice matrix Nash–Sutcliffe Efficiency Total precipitation Liquid precipitation rate Equivalent sphere radius Bulk snow density Dry density of the snowpack Ice density Density of new snow Liquid water density Runoff rate Solid precipitation mass flux Solid precipitation rate Effective saturation degree of the mixture Snow depth Average saturation degree of the mixture Irreducible saturation degree Specific surface area of snow Snow Water Equivalent Total vertical stress of the dry skeleton Air temperature Snow temperature Threshold temperature of melting Mass liquid water content Total volume of the snow mixture Volume of melt–freeze ice Total volume of pores in the snowpack Volume of the porous matrix Volume of the liquid constituent of the mixture Wind speed Elevation

α, α LWC e

η η 0 , f 1 , cη , a η , b η F h hMF hS hW I, I k K KW L LWC, or θ W M MMF MS MW n NSE P p re

ρ ρD ρi ρ NS ρW

R S s S SD Sr Sri SSA SWE

σ

T TS Tτ

W V VMF VP VS VW W z

p is rain precipitation rate (in m/h). We parameterized the melting term following a temperature-index approach [34]. Consequently, M = ρD I[a(T − Tτ )], where I is the product of a first function equal to 1 if T > Tτ , and 0 otherwise, and a second function that tends to 0 with hS [13] (to avoid negative predictions when hS is low), Tτ is the threshold temperature for melting (here equal to 0 °C), and a is a calibration parameter (sometimes referred to as degree–hour parameter in the literature [13,85]). Runoff R is here parametrized following a matrix flow approach [10,72]. This means that

R = ρW · α · KW ,

(4)

where KW is the intrinsic permeability of water in snow, and α is a constant [10]. In particular, α = α  · (5.47 · 105 m−1 s−1 ) [16], where α  is a time conversion constant (e.g. equal to 3600 s/h). Following [10], KW = K · S3 , where K is the intrinsic permeability of snow, in m2 , and S is the effective saturation degree of the mixture, equal to (Sr − Sri )/(1 − Sri ). Sr is the average saturation degree of the porous

4

F. Avanzi et al. / Advances in Water Resources 86 (2015) 1–13

matrix. This is the ratio between liquid water volume and total pores volume. Since the model let the mixture to be over-saturated, we impose Sr = 1 when hW ≥ nhS , and Sr = hW /(nhS ) otherwise. Sri is the irreducible saturation degree. To compute Sri , we consider the formula suggested by [41]: Sri = 0.02ρD /(ρW n). As for K, we adopt here the parameterization proposed by [8]:

K = 3re2 e−0.013·ρD .

(5)

In Eq. (5), re is the equivalent sphere radius, which is defined as 3 2 re = SSA· ρ with respect to snow Specific Surface Area SSA, in m /kg. i

SSA is defined as the area of the ice/air interface in a given snow sample per unit mass [9,42,48,53,69]. This is here estimated adapting the formula proposed by [17]. Accordingly, SSA = −308.2 · ln(ρD ) − 206.0 (SSA is in cm2 g−1 and ρ D in g cm−3 ). For the sake of numerical stability, when Sr > 0.5, Eq. (4) is substituted with a kinematic wave 1.25 . This formula has been proposed by approximation: R = ρW θW hW [13], and it has a marginal impact on snowpack runoff prediction during the snow season. Eq. (4) is constrained by available liquid water and by Sri according to the definition of irreducible water content. Given that MS = ρD hS and MW = ρW hW , Eqs. (2) and (3) can be rewritten, after some algebra, as:

dhS = dt

hS ρD ρNS − I[a(T − Tτ )] s− · ρD ρD dt

dhW = p+ dt

ρD I[a(T − Tτ )] − α · KW ρW

(6)

(7)

We couple these two equations with the momentum balance equation for the dry phase, which is the constituent subjected to deformations (see [13] for details), and with a viscous rheological equation as follows:

d ρD = ρD · dt

σ , η

(8)

where σ is the total vertical stress of the dry skeleton, equal to gρ D hS , g is acceleration due to gravity, and η is viscosity. Here, η = ρ f1 η0 cηD eaη ·(−TS )+bη ρD , following the formulation suggested by [7,77]

for the CROCUS model. With respect to [77], since HyS does not simulate grain shape and/or size, we set to 1 the multiplicative term that relates viscosity with snow microstructure properties (this term is not reported in the definition of η for brevity). Future investigations will address how to include this term explicitly. f1 is a parameter that relates viscosity to liquid water content, while η0 , cη , aη and bη are constants (see [77] for details). TS is the average snow temperature. We predict this quantity following the results by [46]: when T ≥ 0 °C, TS = 0 ◦ C, while when T < 0 °C, a bilinear snow profile is defined (see [13] for details). Since ρ D varies because of new snow events as well, Eq. (8) is modified as follows to account for new events effect on mean dry density, according to its own definition:

d ρD = ρD · dt

σ ρNS − ρD + s. η hS

(9)

See again [13] for the complete derivation. The system of Eqs. (6), (7) and (9) represents the HyS model. It is forced by the hourly input data s, p and T (available at all the stations considered). Depending on the parametrization chosen for fresh snow density, additional data of wind speed may be needed. This formulation needs just one parameter to be calibrated (namely, a). We solve this system numerically [13]. This model is run at hourly time resolution and provides an evaluation of the liquid water content as θW = hW /h. 3.1.2. Melt–freeze temperature-index approach (HyS 2.0 model) Here, we will discuss the inclusion in the HyS model of melt– freeze mass fluxes. This version of the model is indicated as HyS 2.0. Melt–freeze processes are well known to delay liquid water routing in

sub-freezing snow due to the development of ice crusts and/or phase transitions [35,49,50,61,62,76]. Following the observations by [20,79], we will assume that refrozen snow “has a structure and a SSA comparable to that of the parent wet sample, except for effects due to volume expansion” [27]. In particular, we define a generalized dry mass MD = MS + MMF as the sum of the mass of snow grains, i.e. MS (organized in a well structured solid skeleton, and subjected to the dynamics in Eqs. (6) and 9), and the mass of ice produced by the melt–freeze process, i.e. MMF . In this framework, MMF operates as an alternative pore filler, in addition to MW , but with a lower density, hence a greater occupied volume after phase transition. As a first approximation, we assume that this constituent has no effect on SSA dynamics and/or dry density (i.e. the density of snow grains). Accordingly, the dynamics of these variables are predicted using the same equations discussed in Section 3.1.1. MMF occupies a volume equal to VMF , of unitary area and depth equal to hMF (MMF = ρi hMF ). Consequently, h = hS + hMF + hW − nhS  and ρ = (MD + MW )/V = (ρD · hS + ρW · hW + ρi · hMF )/h. The mass balance equations in the integral form for MMF and MW now read

dMMF = +F dt

(10)

dMW =L+M−F dt

(11)

where F is the melt–freeze mass flux. This can be either a positive or negative flux for MMF and MW depending on weather and snow conditions. As M, it is a flux between the two phase. Therefore, it does not imply any mass loss for the mixture as a whole. We model F using a coupled melt–freeze temperature-index approach. If T < 0 °C, we assume F as a mass gain for MMF and a mass loss for MW , provided that MW > 0. If T > 0 °C, we assume F as a mass loss for MMF and a mass gain for MW , provided that MMF > 0. Accordingly, F = −ρW · I · e · a · [T − Tτ ]. 0 ≤ e ≤ 1 is the melt–freeze factor (sometimes also known as degree-day reduction factor for refreezing), and has to be calibrated [65]. e · a has been sometimes referred to as freezing degree hour (or day) factor, depending on model resolution [65]. If T ≥ 0 °C, I = hMF /(hMF + k), while if T < 0 °C I = hW /(hW + k). As for I, this numerical operation is aimed to avoid negative predictions of either hMF and hW . Here, k = 0.01 m. Consequently, HyS 2.0 solves four differential equations. Two are Eqs. (6) and (9), while the modified equations for hW and hMF are as follows:

dhW = p+ dt

ρD I[a(T − Tτ )] + I · e · a · [T − Tτ ] − α · KW ρW

dhMF ρW  I · e · a · [T − Tτ ]. =− dt ρi

(12)

(13)

As for the basic HyS model, HyS 2.0 provides an evaluation of liquid water content at a hourly resolution as θW = hW /h. 3.2. The calibration of the model Both versions of the model need the calibration of a melt degree– hour parameter a at each site. During a given water year, the optimum value of this parameter is selected by using the Least-Squares method between snow depth data and hS . We repeat this operation using each water year one at a time to evaluate parameter variability at the same ¯ location [3]. As a reference value, the median value of a is chosen (a), since this estimate is a resistant measure, unaffected by outliers [47]. If melt–freeze dynamics are considered (HyS 2.0), a second parameter (e) has to be calibrated. For this purpose, we operate here as follows. After calibrating a during a given year (imposing e = 0), we use the Least-Squares method between modeled values of θ W and observation of bulk volumetric LWC in snow pits to estimate e. As for

F. Avanzi et al. / Advances in Water Resources 86 (2015) 1–13

5

a

a

b

c

b Fig. 2. NSE coefficients of snow depth (black symbols), SWE (red symbols) and daily runoff (blue symbols) over the entire dataset in Nagaoka (Fig. 2(a)) and Shinjo (Fig. 2(b)), as a function of the calibration year. Crosses refer to HyS, while dots refer to HyS 2.0. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

a

d Fig. 4. An example of input data and models results in Shinjo for the period 2010 2012. In Fig. 4(a), the black line is air temperature, while the blue horizontal line indicates 0 °C. In Fig. 4(b), the blue and the red lines represent models results in terms of snow depth hS using the HyS and HyS 2.0 approach, respectively. In the same figure, the black line reports data of snow depth, while black dots represent measurements of snow depth in snow pits. In Fig. 4(c) and 4(d) a similar color convention is used to compare models results and data in terms of snow density and SWE, respectively. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

¯ by rea, a median value of e is chosen to be representative at a site (e), peating the calibration procedure using each water year one at a time. This sequential procedure is here considered as a first approximation, but it will be improved in the future. 3.3. The evaluation of the model

b

c

d Fig. 3. An example of input data and models results in Nagaoka for the period 2012– 2014. In Fig. 3(a), the black line is air temperature, while the blue horizontal line indicates 0 °C. In Fig. 3(b), the blue and the red lines (overlapping) represent models results in terms of snow depth hS using the HyS and HyS 2.0 approach, respectively. In the same figure, the black line reports data of snow depth, while black dots represent measurements of snow depth in snow pits. In Fig. 3(c) and (d) a similar color convention is used to compare models results and data in terms of snow density and SWE, respectively. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Models performances are firstly assessed on snow depth, SWE and daily runoff data. The main purpose of this evaluation is the assessment of HyS and HyS 2.0 performances over variables that are usually simulated by snowpack models [5,21,22,72]. A daily resolution for runoff clearly displays its variations at the seasonal scale, but it avoids a more complex analysis on sub-daily variations. The evaluation is performed at Nagaoka and Shinjo, where complete data-sets are available to calibrate both versions of the model. As an index of performances, we chose the well-known Nash–Sutcliffe coefficient NSE [60]. At each site, this coefficient has been calculated by considering the entire simulation horizon and a long-term mean of the three observed variables over the entire study period. The calculation is repeated at each site using both versions of the model, and all calibration solutions available (i.e., one solution per year) to evaluate how performances change as a function of both. Note that, in Shinjo, the NSE coefficient for daily runoff has not been calculated since this type of data is available for just two water years out of seven. As a second step, we evaluate the performances of the model in simulating θ W in Nagaoka and Shinjo, where direct measurements of this quantity are available. We consider as a measure of model performances the absolute difference between observed and modeled volumetric LWC. We call this variable LWC . A set of LWC values can be obtained by picking all the available couples over the entire

6

F. Avanzi et al. / Advances in Water Resources 86 (2015) 1–13

a

or in Shinjo. We then evaluate if models are able to detect dry periods during each season. We assume snow as dry at a given instant if, at that time, snow temperature in the pits was negative at all elevations [73]. As a second evaluation, we consider instants when snow pits data indicated a constant temperature of 0 °C for all elevations (isothermal snow cover). In these conditions, it is likely that θ W > 0. Usually, measurements of LWC are rarer than snow depth data. It follows that calibrating e can be sometimes more challenging than choosing a suitable value for a. Therefore, this preliminary evaluation helps to understand whether the transfer of e from gauged to ungauged sites allows to infer some basic properties of a snow cover (i.e., wet or dry conditions).

b

4. Results

Fig. 5. Two examples of data and models results in terms of daily runoff for the water year 2013. Fig. 5(a) depicts a comparison in Nagaoka, while Fig. 5(b) reports a similar evaluation in Shinjo. The blue and the red lines represent models result using the HyS and HyS 2.0 approach, respectively. The black lines report snow lysimeter data. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

simulation horizon. We therefore investigate how this collection varies as a function of the site, the version of the model and the calibration year considered. Results are summarized using box plots. As a third step, we chose to operate an evaluation of the models in Sapporo, where snow pits data do not include measurements of LWC. We run the models using the median values of a that we obtained locally and the median value of e obtained either in Nagaoka

4.1. Snow depth, density and SWE The calibration of the model returned the following median values of parameters. In Nagaoka, a¯ = 4.94 · 10−4 m/h/°C, while e¯ = 0.2. In Shinjo, a¯ = 5.99 · 10−4 m/h/°C, while e¯ = 0.25. In Sapporo, a¯ = 8.37 · 10−4 m/h/°C. The value of a¯ at a given location does not change either using HyS or HyS 2.0 due to the sequential calibration procedure. We report in Fig. 2 the Nash–Sutcliffe coefficients for snow depth, SWE and daily runoff in Nagaoka and Shinjo (over the entire study period), as a function of the calibration year considered to estimate a and e. An example of models result vs. data in Nagaoka is depicted in Fig. 3. Results are plotted using the a¯ and e¯ values that we obtained locally. Note that these plots are only illustrative, while we

Fig. 6. Box plots of LWC values obtained in Nagaoka (a) and Shinjo (b) during the entire period of simulation, as a function of the calibration year used. As an example, the box plot corresponding to the label 97 in Fig. 6(a) represents the distribution of all LWC obtained when using the calibration solution estimated using water year 1997, and so on. LWC is the absolute difference between observed and simulated θ W on a given day. Blue boxes refer to HyS, while red boxes refer to HyS 2.0. Outliers are represented as crosses. Note that any LWC > 0.045 has not been reported for a matter of readability (see Section 4.2). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

F. Avanzi et al. / Advances in Water Resources 86 (2015) 1–13

7

a

b

c

d

e Fig. 7. Some examples of comparison in terms of θ W between models outputs and data in Nagaoka, during different water years. In all the figures, the blue and the red lines represent models result using the HyS and HyS 2.0 approach, respectively. Black dots are measurements.

refer to Fig. 2 for an exhaustive evaluation of models performance. Just three years are reported (2012–2014), but results are exemplary of the other years, too. Note that HyS and HyS 2.0 results in Figs 3(b), (c) and (d) mostly overlap. Any model output obtained when simultaneous hS is lower than 0.01 m has been set to zero. Fig. 4 depicts some exemplary results at Shinjo. Plots refer to the period 2011–2013, and are representative of the entire time-series. The a¯ and e¯ used are those obtained locally. Any model output obtained when simultaneous hS is lower than 0.01 m has been set to zero. In Fig. 5, two comparisons between daily runoff data and models output are reported. As for previous plots, we used at each site the a¯ and e¯ values that we obtained locally. In Fig. 5(a), an example in Nagaoka is reported, while in Fig. 5(b), an example in Shinjo is depicted. In both the cases, the water year 2013 has been used. 4.2. Liquid water content Fig.s 6(a) and (b) report the box plots of LWC in Nagaoka and Shinjo during the entire period of simulation. Results are given as a

function of the calibration year used. Note that any LWC > 0.045 has not been reported for a matter of graphs readability. However, these values were included to calculate box plot statistics. In Nagaoka, LWC is greater than 0.045 in 10/285 cases when using the 1997 calibration solution (∼ 3.5% of the cases), in 4/285 cases when using the 2001 calibration solution (∼ 1.4%), and in 1 case out of 285 for all other calibration solutions. In Shinjo, LWC is always smaller than 0.045. Figs. 7 and 8 report illustrative comparisons between data and models outputs in terms of θ W at Nagaoka. Fig. 9 depicts time-series of data and models outputs during five out of seven years in Shinjo. In all the cases, we used the a¯ and e¯ values obtained at each site to force the model. We report in Fig. 10 the results of the evaluation of the models in Sapporo. We consider here five out of eight years, but results are consistent over the entire period. In all these plots, any model output obtained when simultaneous hS is lower than 0.01 m has been set to zero. Note that the y axis range has been set to 0–0.1 since this is the typical range of variability of observed LWC (see e.g. [23,24,73]). This is also the range of all observed values and of simulations during most of the season. However, modeled θ W reaches 100 vol% any time the snow cover melts completely (due to the phase

8

F. Avanzi et al. / Advances in Water Resources 86 (2015) 1–13

a

b

c

d

e Fig. 8. Some examples of comparison in terms of θ W between models outputs and data in Nagaoka, during different water years. In all the figures, the blue and the red lines represent models result using the HyS and HyS 2.0 approach, respectively. Black dots are measurements. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

change from ice to liquid water). This behavior is observed at the beginning and at the end of each snow season, only. 5. Discussion 5.1. θ W simulation In Nagaoka, independent from calibration solution or model version, median values of LWC are lower than 1 vol% (Fig. 6(a)). Upper quartiles (i.e., the superior side of box plots) are usually close to 1 vol%. Models outputs are similar, since the differences between median LWC one would obtain using HyS and HyS 2.0 span between 0 and 0.0008 depending on the calibration year considered. The mean absolute difference between HyS and HyS 2.0 predictions when us-

ing local a¯ and e¯ is equal to 0.03 vol% over the entire horizon. Outliers are usually higher than 2.5 % but lower than 4 vol%. Such a range for outliers is reduced if compared with the possible range of variation of θ W in seasonal snow. In fact, this variable can span between 0 and > 15 vol% [24]. However, improving the performances of the models on this side is an important direction of future research. Note that the data-set used here includes two water years (i.e. 1997 and 2003) when data of θ W were usually collected on a daily basis. This represents a valuable information since usually excavating snow pits is very time consuming. Therefore, the temporal resolution of LWC data is usually coarser than daily. In Shinjo, considering melt–freeze dynamics improves the performances of the model in estimating θ W . In fact, while median LWC are greater than 1 vol% when using HyS, this variable is consistently

F. Avanzi et al. / Advances in Water Resources 86 (2015) 1–13

9

a

b

c

d

e Fig. 9. Some examples of comparison in terms of θ W between models outputs and data in Shinjo, during different water years. In all the figures, the blue and the red lines represent models result using the HyS and HyS 2.0 approach, respectively. Black dots are measurements. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

lower than 1 vol% when considering HyS 2.0. Upper quartiles are higher than 1 vol% both for HyS and HyS 2.0, but those by HyS 2.0 are smaller than those by HyS (average difference equal to 0.25 vol%). Focusing on the seasonal dynamics (Fig. 9), HyS 2.0 detects the transition between the dry period during January and February and the wet period during March and April. Overall, the number of LWC outliers that one would obtain by using a given calibration solution are less than in Nagaoka (compare Fig. 4(b) with Fig. 4(a)). However, this may be also due to the fact that more data are available in Nagaoka than in Shinjo. In Sapporo, HyS 2.0 predicts a dry condition of the snow cover in 18 out of 22 observed cases, i.e. 80%. On the contrary, HyS cannot detect dry periods correctly. Isothermal conditions of the snow cover (at 0 °C) are characterized by predictions of θ W > 0. An isothermal snow cover at 0 °C is usually prone to moist or wet conditions, but dry snow can be observed even at 0 °C [24], while θ W can be greater than 0 even if the snow cover is not isothermal at 0 °C. This happens, e.g. due to the onset of superficial percolation in snow at sub-freezing conditions [73]. It follows that this evaluation is just a preliminary

test, since data of θ W are absent. Anyway, it suggests that calibration values of e¯ may be transferred from one site to another in order to get some basic information about the wetness conditions of a snow cover. Note that the calibrated e¯ values in Nagaoka and Shinjo are very similar. This provides useful indications for ungauged applications that will be investigated at other sites in the future. According to [23], the absolute error in measuring water content using dielectric methods spans between 0.2 and 0.9 vol%, depending on dielectric constant, density and their uncertainties (see Eq. (2) and Fig. 1 in [23]). Techel and Pielmeier [73] report 0.5 vol% as a general indication of accuracy for dielectric methods, but they also report that 1 vol% is the expected difference between measurements taken using a Denoth meter and a Snow Fork. These values are close to the median values of LWC that we obtained using HyS 2.0 both in Nagaoka and in Shinjo (< 1 vol%). It follows that HyS 2.0 returns performances in estimating observed θ W that are comparable to the accuracy of the instruments used to observe this variable. Therefore, this approach may represent a practical tool in applications where data availability is restricted. As an example, [57] suggest the ratio between θ W

10

F. Avanzi et al. / Advances in Water Resources 86 (2015) 1–13

a

b

c

d

e Fig. 10. Some examples of comparison in terms of θ W between models outputs and data in Sapporo, during different water years. In all the figures, the blue line denotes model results using the HyS approach, while the green and red lines represent models results using the HyS 2.0 approach, and the e¯ value calibrated in Nagaoka and Shinjo, respectively (lines almost overlap). Black dots represent instants when snow was dry, while black vertical lines indicate instants of isothermal snow at 0 °C. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

and 0.03 (i.e., an estimation of the threshold value of the transition between the pendular and the funicular regimes) as a proxy variable to infer wet snow avalanche activity. Moreover, a coupled melt– freeze temperature-index approach has been sometimes applied in the hydrologic modeling of snow-dominated catchment to replicate the seasonal and/or diurnal variability of streamflow [25,65]. However, performances evaluations at the local scale were still missing, to our knowledge, but provided here. Note that in Nagaoka competitive performances can be obtained considering HyS. All figures show that both models versions return highly fluctuating values of θ W at the beginning and/or at the end of each water year. These fluctuations span between 0 and ∼ 9 vol%, but modeled θ W is sometimes higher than 10 vol% (see e.g. Fig. 7(a), (b), (c) and (e), where the y axis range has been fixed to 0–0.1). The temporal scale of this variation is usually short if compared with the duration of a snow season. This behavior is observed during periods of very shallow snow covers that are undergoing melting. As an example, snow depth never exceeded 0.3 m during the long period of θ W fluctuations

in December 2003–January 2004 in Nagaoka (Fig. 7(e)). At the beginning of the season, these fluctuations usually occur when a single snow event creates a snow cover that melts completely before a new event occurs. In these situations, measuring LWC is very challenging. In fact, no measurement is available to assess if these dynamics are reasonable or not. However, we suggest that models outputs that are highly variable within short time periods should be interpreted with caution. Models outputs are sometimes overestimated. This happens e.g. during 2011 in Nagaoka (Fig. 8(b)). Weather conditions during January and February 2011 were characterized by frequent snowfalls. Moreover, average air temperature in January was equal to −0.33 ◦ C. This is the minimum value of that variable during the considered period (values spanning between −0.33 ◦ C and +2.74 ◦ C). A hypothesis is that the calibrated value of e is not adapt to reproduce θ W during such unusual weather conditions. Note that frequent snow events may cause a great vertical variability in snow density and this may hamper the application of Eq. (1). Anyway, the bias between

F. Avanzi et al. / Advances in Water Resources 86 (2015) 1–13

models and data is temporary, since the agreement increases during late February.

5.2. Snow Depth, density, SWE and runoff When considering SD, NSE for HyS (HyS 2.0) are ≥ 0.85 (≥ 0.85) in Nagaoka, and ≥ 0.90 (≥ 0.90) in Shinjo. In terms of SWE, NSE are ≥ 0.77 (≥ 0.77) in Nagaoka, and ≥ 0.86 (≥ 0.85) in Shinjo. In the hydrological literature, a NSE value close to 1 is usually considered as a proxy of good performances. Fig. 3 shows no significant difference between HyS and HyS 2.0 dynamics in a wet-snow area like Nagaoka. In fact, the average absolute difference between the two model outputs is equal to 0.0014 m, 0.64 kg/m3 and 0.94 mm w.e. in terms of snow depth, density and SWE, respectively. In Shinjo, the average absolute difference between the predictions are equal to 0.004 m, 1.63 kg/m3 and 2.70 mm w.e. for snow depth, density and SWE, respectively. SD and SWE peaks are usually higher considering the melt–freeze process than neglecting it. A higher SD is due to the fact that refreezing reduces θ W , while it increases viscosity. This reduces the settling rate (Eq. (8)). Considering melt–freeze dynamics causes an increase in SWE as well due to mass retention and runoff delay. However, predictions of HyS 2.0 are worse than those by HyS during some years in Shinjo (Fig. 2(b)) or in Nagaoka (Fig. 2(a)). The effect of melt–freeze dynamics on SWE dynamics is clearly low, since the mass contributions of liquid water and refrozen ice are small if compared with the mass of snow grains. This means that the detection of any change in SWE due to this process may be affected by instruments resolution [3]. A clear example is the comparison between data and models during water year 2013 in Shinjo: snow pits data are closer to HyS predictions than HyS 2.0, but snow pillows data are in agreement with the latter one. NSE values for daily runoff are more variable, and usually lower than those obtained for SD and SWE. This is widely expected since snowmelt runoff prediction is still a difficult task even for complex, multi-layer models [80]. Fig. 5(a) shows an example of HyS limitations during March and April 2013 in Nagaoka. In particular, modeled runoff at the beginning of March is greater than observed. This overestimation has probably caused an associated underestimation during the second part of that period. The snow pit profile on February 28 (i.e. at the beginning of this period) reports a relatively high snow depth (183 cm), characterized by the existence of a superficial layer with very fine grains and wet conditions (elevation z between 170 and 147 cm, size ≤ 0.2 mm, and θ W (z) = 5.5%). This overlies a thick basal layer with coarser grains (elevation z between 147 and 0 cm, size between 1 and 2 mm). At the boundary between the layers, an abrupt reduction of θ W (z) is observed (θ W (z) = 0.64% below the boundary). In these textural conditions, capillary barriers often develop. These retain water and delay its undisturbed percolation in the snow medium [12,31,56,78,80,83,84]. This may cause a delay in runoff that cannot be simulated easily by a one-layer model. It is worth noting that, during the same period, observations and simulations of θ W are in agreement (Fig. 8(d)). Similar considerations can be made in Shinjo, too. In that case, an overestimation of water discharge is visible at the beginning of March, while models are not able to reproduce the strong peak in measured runoff in April. It is our opinion that the increase in a¯ with increasing latitude of the site is due to different temporal dynamics of melting. In Nagaoka, snow melting occurs during the entire snow season. On the contrary, in Sapporo winter temperature is usually below 0 °C, and most of melting occurs during spring. Winter conditions are usually characterized by a less amount of direct solar radiation than spring. This could affect the a¯ differences among the sites [34].

11

6. Conclusions We investigated the performances of a simple temperature-index model in reproducing the dynamics of bulk volumetric liquid water content of seasonal snow. Three test sites in Japan, subjected to different climatic conditions, were considered. This analysis showed that, when a coupled melt–freeze temperature-index approach is used, median absolute differences between observed and simulated bulk LWC (LWC ) are smaller than 1 vol% at the two sites where measurements of θ W are available, independent from the year used to calibrate the model. When transferring a calibration parameter (i.e., e) where data to calibrate it are not available (in Sapporo), the model predicts correctly a dry conditions in 18 out of 22 observed cases, i.e. 80%. At the same site, observed isothermal conditions of the snow cover (i.e., measurements of temperature in snow pits returned 0 °C at all elevations sampled) correspond to predictions of θ W that are greater than 0. Overall, this analysis showed that HyS 2.0 returns median performances in estimating observed θ W that are of the same order of magnitude of the accuracy of the instruments used to observe this variable (that is, median absolute difference between data and predictions lower than 1 vol%). This proves that a coupled melt– freeze approach may be used to infer an average value of LWC in seasonal snow. However, it is also found that uncalibrated weather or snow conditions (e.g., a heavy snowfall period and/or a marked vertical variability in liquid water content) may cause errors in the prediction. Although the impact of these errors is usually restricted in time, additional tests will be run in the future to improve model performances.

Acknowledgments We would like to thank the staff of the Snow and Ice Research Center, National Research Institute for Earth Science and Disaster Prevention, for helpful discussions. FA is grateful for the support received during his research period at the Snow and Ice Research Center. Gratitude is due to Masaki Nemoto and Kenji Kosugi (Shinjo Cryospheric Environment Laboratory, Snow and Ice Research Center, NIED) for Shinjo observations, and to Teruo Aoki (Meteorological Research Institute, Tsukuba-shi, Ibaraki-ken) and Sumito Matoba (Institute of Low Temperature Science, Hokkaido University) for Sapporo observations. Data in Sapporo were obtained with the support by the Grant for Joint Research Program, Institute of Low Temperature Science, Hokkaido University, and the Japan Society for the Promotion of Science (JSPS), Grant-in-Aid for Scientific Research (S), no. 23221004. We would like to thank the Editor, Prof. Paolo D’Odorico, and two anonymous reviewers, for valuable comments on the manuscript.

References [1] Aoki T, Motoyoshi H, Kodama Y, Yasunari TJ, Sugiura K. Variations of the snow physical parameters and their effects on albedo in Sapporo, Japan. Ann Glaciol 2007;46:375–81. http://dx.doi.org/10.3189/172756407782871747. [2] Avanzi F, Caruso M, Jommi C, De Michele C, Ghezzi A. Continuous-time monitoring of liquid water content in snowpacks using capacitance probes: a preliminary feasibility study. Adv Water Resour 2014;68:32–41. http://dx.doi.org/10.1016/j. advwatres.2014.02.012. [3] Avanzi F, De Michele C, Ghezzi A, Jommi C, Pepe M. A processing-modeling routine to use snotel hourly data in snowpack dynamic models. Adv Water Resour 2014;73:16–29. http://dx.doi.org/10.1016/j.advwatres.2014.06.011. [4] Baggi S, Schweizer J. Characteristics of wet-snow avalanche activity: 20 years of observations from a high alpine valley (Dischma, Switzerland). Nat Hazards 2008;50:97–108. http://dx.doi.org/10.1007/s11069-008-9322-7. [5] Bartelt P, Lehning M. A physical snowpack model for the Swiss avalanche warning part i: numerical model. Cold Regions Sci Technol 2002;35:123–45. http://dx.doi. org/10.1016/S0165-232X(02)00074-5. [6] Brun E. Investigation on wet-snow metamorphism in respect of liquid-water content. Ann Glaciol 1989;13:22–6.

12

F. Avanzi et al. / Advances in Water Resources 86 (2015) 1–13

[7] Brun E, David P, Sudul M, Brunot G. A numerical model to simulate snow-cover stratigraphy for operational avalanche forecasting. J Glaciol 1992;38(128):13–22. [8] Calonne N, Geindreau C, Flin F, Morin S, Lesaffre B, du Roscoat SR, et al. 3-d image-based numerical computations of snow permeability: links to specific surface area, density, and microstructural anisotropy. Cryosphere 2012;6(5):939–51. http://dx.doi.org/10.5194/tc-6-939-2012. [9] Carmagnola CM, Morin S, Lafaysse M, Domine F, Lesaffre B, Lejeune Y, et al. Implementation and evaluation of prognostic representations of the optical diameter of snow in the surfex/isba-crocus detailed snowpack model. Cryosphere 2014;8(2):417–37. http://dx.doi.org/10.5194/tc-8-417-2014. [10] Colbeck SC. A theory of water percolation in snow. J Glaciol 1972;11(63):369–85. [11] Colbeck SC. Theory of metamorphism of wet snow. Technical report. Hanover, NH, USA: Cold Regions Research and Engineering Laboratory; 1973. [12] Coleou C, Lesaffre B. Irreducible water saturation in snow: experimental results in a cold laboratory. Ann Glaciol 1998;26:64–8. [13] De Michele C, Avanzi F, Ghezzi A, Jommi C. Investigating the dynamics of bulk snow density in dry and wet conditions using a one-dimensional model. Cryosphere 2013;7(2):433–44. http://dx.doi.org/10.5194/tc-7-433-2013. [14] Denoth A. The pendular-funicular liquid transition and snow metamorphism. J Glaciol 1980;28(99):357–64. [15] Denoth A. An electronical device for long-term snow wetness recording. Ann Glaciol 1994;19:104–6. [16] DeWalle DR, Rango A. Principles of snow hydrology. Cambridge University Press; 2011. [17] Domine F, Taillandier A-S, Simpson WR. A parameterization of the specific surface area of seasonal snow for field use and for models of snowpack evolution. J Geophys Res: Earth Surf. 2007;112 (F2)(ISSN 2156-2202). http://dx.doi.org/10.1029/ 2006JF000512. [18] Domine F, Morin S, Brun E, Lafaysse M, Carmagnola CM. Seasonal evolution of snow permeability under equi-temperature and temperature-gradient conditions. Cryosphere 2013;7(6):1915–29. http://dx.doi.org/10.5194/tc-7-1915-2013. [19] Eiriksson D, Whitson M, Luce CH, Marshall HP, Bradford J, Benner SG, et al. An evaluation of the hydrologic relevance of lateral flow in snow at hillslope and catchment scales. Hydrol Proc 2013;27:640–54. http://dx.doi.org/10.1002/hyp. 9666. [20] Erbe EF, Rango A, Foster J, Josberger EG, Pooley C, Wergey WP. Collecting, shipping, storing, and imaging snow crystals and ice grains with low-temperature scanning electron microscopy. Microsc Res Tech 2003;62:19–32. http://dx.doi. org/10.1002/jemt.10383. [21] Essery R, Martin E, Douville H, Fernández A, Brun E. A comparison of four snow models using observations from an alpine site. Clim Dyn 1999;15(8):583–93. http://dx.doi.org/10.1007/s003820050302. [22] Essery R, Morin S, Lejeune Y, Ménard CB. A comparison of 1701 snow models using observations from an alpine site. Adv Water Resour 2013;55:131–48. http://dx. doi.org/10.1016/j.advwatres.2012.07.013. [23] Fierz C, Föhn MB. Long-term observation of the water content of an alpine snowpack. In: Proceedings of the International Snow Science Workshop 1994; 1994. Snowbird, Utah, USA. [24] Fierz C, Armstrong RL, Durand Y, Etchevers P, Greene E, McClung DM, et al. The international classification for seasonal snow on the ground Technical report. UNESCO-IHP, Paris: IHP-VII Technical Documents in Hydrology N83, IACS Contribution N1; 2009. [25] Formetta G, Kampf SK, David O, Rigon R. Snow water equivalent modeling components in newage-jgrass. Geosci Model Dev. 2014;7(3):725–36. http://dx.doi.org/ 10.5194/gmd-7-725-2014. [26] Gabbi J, Huss M, Bauder A, Cao F, Schwikowski M. The impact of saharan dust and black carbon on albedo and long-term mass balance of an alpine glacier. Cryosphere 2015;9(4):1385–400. http://dx.doi.org/10.5194/tc-9-1385-2015. [27] Gallet J-C, Domine F, Dumont M. Measuring the specific surface area of wet snow using 1310 nm reflectance. Cryosphere 2014;8(4):1139–48. http://dx.doi.org/10. 5194/tc-8-1139-2014. [28] Heilig A, Mitterer C, Schmid L, Wever N, Schweizer J, Marshall H-P, Eisen O. Seasonal and diurnal cycles of liquid water in snow - Measurements and modeling. J Geophysical Research Earth Surface 2015. http://dx.doi.org/10.1002/ 2015JF003593. [29] Heilig A, Eisen O, Schneebeli M. Temporal observations of a seasonal snowpack using upward-looking GPR. Hydrol Process 2010;24:3133–45. http://dx.doi.org/ 10.1002/hyp.7749. [30] Hirashima H, Nishimura K, Baba E, Hachikubo A, Lehning M. Snowpack model simulations for snow in Hokkaido, Japan. Ann Glaciol 2004;38:123–9. http://dx. doi.org/10.3189/172756404781815121. [31] Hirashima H, Yamaguchi S, Sato A, Lehning M. Numerical modeling of liquid water movement through layered snow based on new measurements of the water retention curve. Cold Regions Sci Technol 2010;64:94–103. http://dx.doi.org/10. 1016/j.coldregions.2010.09.003. [32] Hirashima H, Yamaguchi S, Katsushima T. A multi-dimensional water transport model to reproduce preferential flow in the snowpack. Cold Regions Sci Technol 2014;108:80–90. http://dx.doi.org/10.1016/j.coldregions.2014.09.004. [33] Hirashima H, Yamaguchi S, Kosugi K, Nemoto M, Aoki T, Matoba S. Validation of the snowpack model using snow pit observation data (in Japanese with english abstract). Seppyo 2015;77(1):5–16. [34] Hock R. Temperature index melt modelling in mountain areas. J Hydrol 2003;282:104–15. http://dx.doi.org/10.1016/S0022-1694(03)00257-9.

[35] Illangasekare TH, Walter RJ Jr, Meier MF, Pfeffer WT. Modeling of meltwater infiltration in subfreezing snow. Water Resour Res 1990;26(5):1001–12. http://dx. doi.org/10.1029/WR026i005p01001. [36] Ishizaka M. New categories for the climatic division of snowy areas in Japan. Ann Glaciol 1998;26:131–7. [37] Kaempfer TU, Schneebeli M. Observation of isothermal metamorphism of new snow and interpretation as a sintering process. J Geophys Res 2007;112:D24101. http://dx.doi.org/10.1029/2007JD009047. [38] Kajikawa M. Relationship between new snow density and shape of snow crystals (in japanese). Seppyo 1989;51(3):178–83. [39] Katsushima T, Kumakura T, Takeuchi Y. A multiple snow layer model including a parameterization of vertical water channel process in snowpack. Cold Regions Sci Technol 2009;59:143–51. http://dx.doi.org/10.1016/j.coldregions.2009.09.002. [40] Katsushima T, Yamaguchi S, Kumakura T, Sato A. Experimental analysis of preferential flow in dry snowpack. Cold Regions Sci Technol 2013;85:206–16. http://dx. doi.org/10.1016/j.coldregions.2012.09.012. [41] Kelleners TJ, Chandler DG, McNamara JP, Gribb MM, Seyfried MS. Modeling the water and energy balance of vegetated areas with snow accumulation. Vadose Zone J 2009;8(4):1013–30. http://dx.doi.org/10.2136/vzj2008.0183. [42] Kerbrat M, Pinzer B, Huthwelker T, Gäggeler HW, Ammann M, Schneebeli M. Measuring the specific surface area of snow with x-ray tomography and gas adsorption: comparison and implications for surface smoothness. Atmos Chem Phys 2008;8(5):1261–75. http://dx.doi.org/10.5194/acp-8-1261-2008. [43] Kienzle SW. A new temperature based method to separate rain and snow. Hydrol Proc 2008;22:5067–85. http://dx.doi.org/10.1002/hyp.7131. [44] Kirchner HOK, Michot G, Narita H, Suzuki T. Snow as a foam of ice: plasticity, fracture and the brittle-to-ductile transition. Philos Mag A 2001;81(9):2161–81. http://dx.doi.org/10.1080/01418610108217141. [45] Koch F, Prasch M, Schmid L, Schweizer J, Mauser W. Measuring snow liquid water content with low-cost GPS receivers. Sensors 2014;14:20975–99. http://dx.doi. org/10.3390/s141120975. [46] Kondo J, Yamazaki T. A prediction model for snowmelt, snow surface temperature and freezing depth using a heat balance method. J Appl Meteorol 1990;29:375– 84. http://dx.doi.org/10.1175/1520-0450(1990)0290375:APMFSS2.0.CO;2. [47] Kottegoda NT, Rosso R. Applied statistics for civil and environmental engineers. McGraw-Hill; 2008. [48] Legagneux L, Cabanes A, Dominé F. Measurement of the specific surface area of 176 snow samples using methane adsorption at 77 K. J Geophys Res: Atmospheres 2002;107(D17). ACH 5–1–ACH 5–15, ISSN 2156-2202, http://dx.doi.org/10.1029/ 2001JD001016. [49] Marsh P, Woo M-K. Wetting front advance and freezing of meltwater within a snow cover 1. observations in the canadian arctic. Water Resour Res 1984;20(12):1853–64. http://dx.doi.org/10.1029/WR020i012p01853. [50] Marsh P, Woo M-K. Wetting front advance and freezing of meltwater within a snow cover 2. a simulation model. Water Resour Res 1984;20(12):1865–74. http://dx.doi.org/10.1029/WR020i012p01865. [51] Marsh P, Woo MK. Meltwater movement in natural heterogeneous snow covers. Water Resour Res 1985;21(11):1710–16. http://dx.doi.org/10.1029/ WR021i011p01710. [52] Marshall HP, Conway H, Rasmussen LA. Snow densification during rain. Cold Regions Sci Technol 1999;30:35–41. http://dx.doi.org/10.1016/S0165-232X(99) 00011-7. [53] Matzl M, Schneebeli M. Measuring specific surface area of snow by nearinfrared photography. J Glaciol 2006;52(179):558–64. http://dx.doi.org/10.3189/ 172756506781828412. [54] Mitterer C, Schweizer J. Analysis of the snow-atmosphere energy balance during wet-snow instabilities and implications for avalanche prediction. Cryosphere 2013;7(1):205–16. http://dx.doi.org/10.5194/tc-7-205-2013. [55] Mitterer C, Heilig A, Schweizer J, Eisen O. Upward-looking ground-penetrating radar for measuring wet-snow properties. Cold Regions Sci Technol 2011;69:129– 38. http://dx.doi.org/10.1016/j.coldregions.2011.06.003. [56] Mitterer C, Hirashima H, Schweizer J. Wet-snow instabilities: comparison of measured and modelled liquid water content and snow stratigraphy. Ann Glaciol 2011;52(58):201–8. http://dx.doi.org/10.3189/172756411797252077. [57] Mitterer C, Techel F, Fierz C, Schweizer J. An operational supporting tool for assessing wet-snow avalanche danger. In: International Snow Science Workshop Grenoble Chamonix Mont-Blanc - 2013; 2013. [58] Morland LW, Kelly RJ, Morris EM. A mixture theory for a phase-changing snowpack. Cold Regions Sci Technol 1990;17:217–85. [59] Nakamura K, Mochizuki S, Kosugi K, Nemoto M, Sato K, Abe O. Meteorological, snowfall and snow cover data observed at shinjo (2013/14 winter). Technical report. Technical Note of the National Research Institute for Earth Science and Disaster Prevention; 2015.390, p. 47 (in Japanese with English abstract). [60] Nash JE, Sutcliffe JV. River flow forecasting through conceptual models part i - a discussion of principles. J Hydrol 1970;10:282–90. [61] Pfeffer WT, Humphrey NF. Determination of timing and location of water movement and ice-layer formation by temperature measurements in sub-freezing snow. J Glaciol 1996;42(141):292–304. [62] Pfeffer WT, Illangasekare TH, Meier MF. Analysis and modeling of melt-water refreezing in dry snow. J Glaciol 1990;36(123):238–46. [63] Pinzer BR, Schneebeli M. Snow metamorphism under alternating temperature gradients: Morphology and recrystallization in surface snow. Geophys Res Lett 2009;36:L23503. http://dx.doi.org/10.1029/2009GL039618.

F. Avanzi et al. / Advances in Water Resources 86 (2015) 1–13 [64] Saito K, Yamaguchi S, Iwata H, Harazono Y, Kosugi K, Lehning M, et al. Climatic physical snowpack properties for large-scale modeling examined by observations and a physical model. Polar Sci 2012;6:79–95. http://dx.doi.org/10.1016/j.polar. 2012.02.003. [65] Schaefli B, Nicótina L, Imfeld C, Ronco PD, Bertuzzo E, Rinaldo A. Sehrecho v1.0: a spatially explicit hydrologic response model for ecohydrologic applications. Geosci Model Dev. 2014;7(6):2733–46. http://dx.doi.org/10.5194/ gmd-7-2733-2014. [66] Schleef S, Löwe H, Schneebeli M. Hot-pressure sintering of low-density snow analyzed by X-ray microtomography and in situ microcompression. Acta Mater 2014;71:185–94. http://dx.doi.org/10.1016/j.actamat.2014.03.004. [67] Schmid L, Koch F, Heilig A, Prasch M, Eisen O, Mauser W, et al. A novel sensor combination (upGPR-GPS) to continuously and nondestructively derive snow cover properties. Geophys Res Lett 2015;42:1–9. http://dx.doi.org/10.1002/ 2015GL063732. [68] Schneebeli M. Development and stability of preferential flow paths in a layered snowpack. In: Blogeochemistry of Seasonally Snow-Covered Catchments (Proceedings of a Boulder Symposium); 1995.Number IAHS Publ. no. 228. [69] Schneebeli M, Sokratov SA. Tomography of temperature gradient metamorphism of snow and associated changes in heat conductivity. Hydrol Process 2004;18:3655–65. http://dx.doi.org/10.1002/hyp.5800. [70] Shimizu H. Air permeability of deposited snow. Technical report. Contributions from the Institute of Low Temperature Science; 1970. A22, 1–32. [71] Sturm M, Holmgren J, Liston GE. A seasonal snow cover classification system for local to global applications. J Climate 1995;8:1261–83. http://dx.doi.org/10.1175/ 1520-0442(1995)0081261:ASSCCS2.0.CO;2. [72] Tarboton DG, Luce CH. Utah energy balance snow accumulation and melt model (UEB), computer model technical description and users guide. Technical report. Utah Water Research Laboratory Utah State University and USDA Forest Service; 1996. [73] Techel F, Pielmeier C. Point observations of liquid water content in wet snow investigating methodical, spatial and temporal aspects. Cryosphere 2011;5(2):405–18. http://dx.doi.org/10.5194/tc-5-405-2011. [74] Techel F, Pielmeier C, Schneebeli M. Microstructural resistance of snow following first wetting. Cold Regions Sci Technol 2011;65:382–91. http://dx.doi.org/10. 1016/j.coldregions.2010.12.006.

13

[75] Tedesco M, Kim EJ, England AW, De Roo RD, Hardy JP. Brightness temperatures of snow melting/refreezing cycles: observations and modeling using a multilayer dense medium theory-based model. IEEE Trans Geosci Remote Sens 2006;44(12):3563–73. http://dx.doi.org/10.1109/TGRS.2006.881759. [76] Tseng PH, Illangasekare TH, Meier MF. Modeling of snow melting and uniform wetting front migration in a layered subfreezg snowpack. Water Resour Res 1994;30(8):2363–76. http://dx.doi.org/10.1029/94WR00764. [77] Vionnet V, Brun E, Morin S, Boone A, Faroux S, Moigne PL, et al. The detailed snowpack scheme crocus and its implementation in surfex v7.2. Geosci Model Dev 2012;5(3):773–91. http://dx.doi.org/10.5194/gmd-5-773-2012. [78] Waldner PA, Schneebeli M, Schultze-Zimmermann U, Flühler H. Effect of snow structure on water flow and solute transport. Hydrol Process 2004;18(7):1271– 90. http://dx.doi.org/10.1002/hyp.1401. [79] Wergin WP, Rango A, Erbe EF. Observations of snow crystals using lowtemperature scanning electron microscopy. Scanning 1995;17:41–9. [80] Wever N, Fierz C, Mitterer C, Hirashima H, Lehning M. Solving Richards equation for snow improves snowpack meltwater runoff estimations in detailed multilayer snowpack model. Cryosphere 2014;8(1):257–74. http://dx.doi.org/10.5194/ tc-8-257-2014. [81] Yamaguchi S, Sato A, Lehning M. Application of the numerical snowpack model (snowpack) to the wet-snow region in Japan. Ann Glaciol 2004;38:266–72. [82] Yamaguchi S, Nakai S, Iwamoto K, Sato A. Influence of anomalous warmer winter on statistics of measured winter precipitation data. J Appl Meteorol Clim 2009;48:2403–9. http://dx.doi.org/10.1175/2009JAMC2008.1. [83] Yamaguchi S, Katsushima T, Sato A, Kumakura T. Water retention curve of snow with different grain sizes. Cold Regions Sci Technol 2010;64:87–93. http://dx.doi. org/10.1016/j.coldregions.2010.05.008. [84] Yamaguchi S, Watanabe K, Katsushima T, Sato A, Kumakura T. Dependence of the water retention curve of snow on snow characteristics. Ann Glaciol 2012;53(61):6–12. http://dx.doi.org/10.3189/2012AoG61A001. [85] Yamaguchi S, Hirashima H, Sato A. Relationship between preferential flow and conditions of water input in snow cover. In: International Snow Science Workshop Grenoble Chamonix Mont-Blanc - 2013; 2013. [86] Zermatten E, Schneebeli M, Arakawa H, Steinfeld A. Tomography–based determination of porosity, specific area and permeability of snow and comparison with measurements. Cold Regions Sci Technol 2014;97:33–40. http://dx.doi.org/ 10.1016/j.coldregions.2013.09.013.