Snow water equivalent modeling components in NewAge-JGrass

The paper presents a package of a modified temperature index based snow water equivalent model as part of the hydrological modeling system NewAge-JGrass. Three temperature-based snow models are integrated in the NewAge-JGrass modeling system and use many of its com5 ponents such as those for radiation balance (SWRB), kriging (KRIGING), automatic calibration algorithms (particle swarm optimization), and tests of goodness of fit (NewAgeV), to build suitable modelling solutions (MS). Similarly to all the NewAge-JGrass components, the models can be exe10 cuted both in raster and in vector mode. The simulation time step can be daily, hourly or sub-hourly, depending on user needs and availability of input data. The MS are applied on the Cache la Poudre river basin (CO, USA) using three test applications. First, daily snow water equivalent is simulated 15 for three different measurement stations for two snow model formulations. Second, hourly snow water equivalent is simulated using all the three different snow model formulations. Finally a raster mode application is performed to compute snow water equivalent maps for the whole Cache la Poudre 20 basin. In all the applications the model performance is satisfactory in terms of goodness of fit relative to measured snow water equivalent time series and the results, and the differences in performances of the different modelling solutions are discussed. 25


Introduction
The process-based distributed approach is the most complete method of simulating snowpack evolution.This solution has reached maturity and was pursued successfully with many recent models including CROCUS (Brun et al. (1992)), Alpine3D (Lehning et al. (2006)), GEOtop (Rigon et al. (2006); Endrizzi et al. (2013); Endrizzi (2007); Dall'Amico  2011)), ISNOBAL (Marks et al. (1999)), and UEB (Tarboton et al. (1996)).These models simulate snow accumulation and ablation using the energy budget and may also include ancillary modeling of blowing snow and other features required to reproduce the full set of thermodynamic 40 quantities representative of snowpack state.However, realtime modeling with data assimilation and parameter calibration may require that a forecasting simulation be generated in few minutes, and this can be accomplished only with simpler models.Simpler models may be limited to forecast-45 ing just the snow water equivalent (SWE, the mass of liquid water in the snowpack) and not other variables such as snow depth and density.One early example of a simple snow accumulation and ablation model is the Snowmelt Runoff Model (SRM) by Martinec (1975).This model was applied 50 to hundreds of basins with reasonable success Martinec et al. (1983) and Martinec et al. (1994).SRM is a linear model in which the independent variables are average daily temperature and an estimate of the catchment area covered by snow.The snow-covered area can be determined from airborne or 55 satellite remote sensing data, and loss of snow cover was then simulated based on a temperature index.Simulations were typically run at a daily time step.
In other studies (Cazorzi and Dalla Fontana (1996); Hock (1999)) the temperature index snow modeling was improved 60 by incorporating a radiation term in addition to temperature.In Cazorzi and Dalla Fontana (1996) the radiation term is an energy index computed for each pixel of the grid as shortwave solar radiation integrated over time, as explained in section 2. In Hock (1999) the melt factor depends on the value 65 of the clear sky solar radiation, following on studies by Kustas et al. (1994) and Brubaker et al. (1996).Hock's model Formetta et al.: Snow water equivalent modeling component in NewAge-JGrass depends on two separate terms: a constant value (melt coefficient) and a value function of the potential solar radiation (radiation coefficient).A third temperature-based snow modeling approach was presented by Tobin et al. (2012) who proposed to use a varying degree-day factor throughout the day to improve simulation of snowmelt rates at sub-daily time steps as a component of a runoff model.
In this paper we implement three of these temperaturebased snow models: a temperature index (C1), Cazorzi and Dalla Fontana's model (C2) and Hock's model (C3) of snow water equivalent, that estimates SWE from spatially distributed radiation and temperature.They are provided as Object Modeling System version 3 (OMS3, David et al. (2013)) components and integrated with the other components of the JGrass-NewAGE system.This system provides an optimal framework for comparing modelling solutions, as all the ancillary tools used remain unchanged when switching from one SWE model to the other.The model components can then be executed using OMS3 implicit parallelism to improve computational efficiency in multicore or multiprocessor machines.The paper is organized as follows: section 2 presents the models' equations; section 3 contains a general description of the JGrass-NewAge system, and section 4 contains a test of the model for an example basin.

The NewAge-SWE Component
The snow water equivaltent modeling components in NewAge are built following the conceptual scheme presented in Kokkonen et al. (2006), varying the contents of the snowpack mass balance equation.In particular the new additions are: -Snow melt is simulated with three different temperature-based solutions; -The separation of rain from snow precipitation uses a smoothing function based on air temperature rather than a threshold air temperature.This approach addresses problems found in prior research by Kavetski et al. (2006) who found that threshold temperatures for precipitation could generate extremely non-smooth parameter surfaces during automatic calibration procedures.
In the next subsection the main algorithms of the model are described with more detail.

Mass Balance
The snowpack mass balance is computed as follows.For the water equivalent of ice (M i [L]): and for liquid water (M w [L]) in the snowpack: Eq.( 1) represents time-varying ice in the snowpack as the 115 sum of snowfall, P s , and freezing water, F, minus melt, M (all expressed as snow water equivalent).Eq.( 2) represents timevarying liquid water in the snowpack as the sum of the rainfall, P r , and melt water minus freezing water.If liquid water M w exceeds liquid water-retention capacity of the snowpack 120 (M max [mm]), the surplus becomes snowmelt discharge q m [L/T].The liquid water retention capacity of a snowpack is related to the ice content by a linear relationship depending on the coefficient α l [-], eq.( 3) 125 Kokkonen et al. (2006) computed these mass balance equations at a daily time step, but here the time step can vary depending on the time resolution of input data.

Type of precipitation
The first hydrological process simulated is the discrimina-130 tion between rainfall and snowfall considering that the two forms of precipitation appear as distinct in equations ( 2) and (1).Usually only precipitation totals and air temperature are available from meteorological stations.A common procedure for separating rain and snow is to use a threshold air temper-135 ature T s : all the precipitation is considered snow if the air temperature for the time interval is less than or equal to T s ; all the precipitation is considered to be rain if air temperature is greater than T s .As proposed in Kavetski et al. (2006) to avoid problems for parameter calibration, a smoother filter 140 for thresholds is applied, and the algorithm to discriminate between rainfall and snowfall can be described as follows: where: P [L/T] is measured precipitation, P r [L/T] is the rainfall precipitation, P s [L/T] is the snowfall precipitation, T s 145 [C] is the threshold temperature and m 1 [-] is the parameter controls the degree of smoothing (if m 1 → 0 threshold behavior is simulated).The two coefficients α r , and α s adjust for measurement errors for rain and snow.Because different values for different climate regions have been found 150 in prior studies (Forland et al. (1996); Rubel and Hantel (1999); Michelson (2004)), in the model the two coefficients are considered parameters and are therefore calibrated.

Snow melt fluxes
The model includes three snow melt formulations.The user 155 is able to select the any of these depending on the site characteristics and data availability.The first melt component (C1) is a traditional temperature index method where the snow melt rates depends only on air temperature: is the melt factor, and T [C] is the air temperature.The model can be used either at hourly or daily time steps, if the parameters are calibrated accordingly.
The second snow melt component (C2) is based on the approach presented in Cazorzi and Dalla Fontana (1996): the melt rate is a function of both shortwave radiation and air temperature.The equation for melt during the day is The equation for the melt process during the night is: where: ] is an energy index, and V s [-] is the sky view factor.This energy index is the mean energy from shortwave radiation over a given period at a certain point, and can be variable in space.In practice, the shortwave direct and diffuse solar radiation is estimated by means of an appropriate tool.In this paper, as presented in Cazorzi and Dalla Fontana (1996) five EI maps are computed starting from 21 December (winter solstice) to the middle of February, March, April, May and June.Different time intervals could be selected depending on user needs.During the night the snow melt is a function of the energy index minimum value of the considered map, as presented in Cazorzi and Dalla Fontana (1996).
The third snow melt component (C3) is based on the formulation presented in Hock (1999).Unlike C2, where the energy index is variable in space but integrated over time, C3 requires the computation of the solar energy for each timestep of the simulation.The melt formulation is: where R s [E] is the incoming (beam plus diffuse) solar radiation received by the pixel and computed using the model presented in Formetta et al. (2013b) The shortwave radiation model from Formetta et al. (2013b), unlike the radiation model presented in the original Hock's formulation, is able to account for shadow effects, complex topography, and compute diffuse radiation.

Freezing
The rate of freezing F used in the mass balance equations is a linear function of air temperature when the air temperature is less then the melting temperature T m , as presented in eq.( 9) where F [L T −1 ] is the freezing rate, and α f [L C −1 t −1 ] is 205 the freezing factor.
For both melt and freezing, if the model is run at a daily time step temperature is the mean daily temperature.If it is applied at and hourly time step, temperature is the mean hourly temperature.The JGrass-NewAge system (Formetta ( 2013)) provides a pool of model components that can be connected and exchanged at run-time.A working set of components constitutes a modelling solution (MS) which is usually set up for a 215 particular purpose or set of simulations.A MS is actually tied together by means of a scripting language (a domain specific language or DSL), and the scripts can be stored together with the input data to preserve memory of a certain simulation set, which can then be easily reproduced and inspected by third 220 parties.JGrass-NewAGE includes components simulate various hydrological processes such as: the space-time structure of precipitation, (KRIGING) shortwave and long wave radiation balance, (SWRB and LWRB), Formetta et al. (2013b).
aggregation and routing of flows in channel (Routing), Formetta et al. (2011) 230 Finally, it also includes three different automatic calibration algorithms: -Particle Swarm Optimization component (Eberhart and Shi (2001))

235
-LUCA (Let Us CAlibrate) component (Hay et al. (2006)) -DREAM component (Vrugt et al. (2009)) The system is based on a hillslope-link geometrical partition of the landscape, so the basic unit for the water budget eval-240 uation is the hillslope.Each hillslope, rather than a cell or a pixel, drains into a single associated link.The model requires interpolation of the meteorological forcing data (air temperature, precipitation, relative humidity) for each hillslope.This operation can be handled by a deterministic inverse distance weighted algorithm (Cressie (1992); Lloyd (2005)), kriging (Goovaerts (1997)) or detrended kriging as in Garen et al. (1994) and Garen and Marks (2005).The radiation model (Formetta et al. (2013b)) implements algorithms that take into account shadows and complex topography.Shortwave radiation under generic sky conditions (all-sky) is computed according to Helbig et al. ( 2010) and using different parameterizations choices such as Erbs et al. (1982), Reindl et al. (1990) and Orgill and Hollands (1977).The longwave radiation budget is based on Brutsaert (1982) and Brutsaert (2005).
All modeling components (including those not described here) can be calibrated using one of the automatic calibration algorithms implemented: the Particle Swarm Optimization algorithm, LUCA and DREAM.Verification of each model component's behavior is eventually tested with the use of NewAge-V (Verification), which provides some of the classical indices of goodness of fit such as: Nash-Sutcliffe, Percent bias, Index of agreement and Kling Gupta efficiency, all defined in Appendix A. The complete interoperating set of components available so far can be seen in fig.(11).
The snow melt model components, SWE-C, are perfectly integrated in the NewAge System as presented in fig.(12).It uses the kriging tools for spatial interpolation of temperature and precipitation and another interpolation method, JAMI, presented in Formetta (2013) for temperature interpolation.Like the interpolation algorithms, SWE-C can be applied both to raster grids and for individual points.SWE-C also uses the NewAge short wave radiation component to estimate the maps of accumulated energy in different periods of the year based on topography, shadow, and cloud cover.The SWE-C outputs could be raster maps or time-series of snow water equivalent and snow melt for any point within the domain.If coupled with runoff modeling, these points could be centroids of hillslopes.The SWE-C component could be connected to the NewAge and OMS3 calibration algorithm to estimate the best model parameters values.
The MS shown in fig.( 12) can be further connected to other available components to obtain an estimation of the runoff, although demonstration of this application is not the goal of this paper.

Sites and data description
To test the performance of SWE-C, the model is applied in the upper Cache la Poudre River basin, located in the Rocky Mountains of northern Colorado and southern Wyoming, USA.This basin is 2700 km 2 has elevations ranging from 1590-4125 m, with mean annual precipitation ranging from 330mm at lower elevations to 1350mm at the highest elevations (Richer et al. (2013) 12) using the Particle Swarm calibration algorithm.The first year of the available time series was selected as calibration period for each station.The "optimal" parameter set estimated in the calibration period was used for the model application in the remaining part of the available 325 time series (verification period).The Kling-Gupta Efficiency (KGE), see Appendix, presented in Gupta et al. (2009) was selected as the calibration objective function.The appendix also describes the motivation for using the goodness of fit indices used in the paper and presents equations for each: 330 Nash-Sutcliffe Efficiency (NSE), Percentage Bias (PBIAS) and Index of Agreement (IOA).

Test 1: Model calibration and verification at daily time-step.
In this application models C1 and C2 were calibrated and 335 verified.Model C3 was not applied in this case because the SWRB component needed air temperature and relative humidity at an hourly time step to compute the incoming solar radiation.As specified in Formetta et al. (2013b), this is necessary for modeling the current atmospheric transmittances.

340
This information was not available at every station.Tables ( 12) and ( 13) report the optimal parameter sets of the models C1 and C2 respectively for each of the three locations where they were applied.Tables ( 14) and ( 15) show for the calibration period and 345 for entire simulation period the goodness of fit values for the models C1 and C2 respectively and for the three SNOTEL stations.Model C2 performs better than the classical temperature index model both in calibration and in verification period in each of the three locations.
Values of the objective function, KGE, were around 0.90 for the model C1 and around 0.95 for model C2, for the calibration periods.In verification period values decreased to 0.80 for C1 model and 0.90 for C2.Similarly, values of the other performance metrics declined from the calibration to verification time periods, but they all are within the range specified for "good" model performance according the guidelines presented in Stehr et al. (2008) and Van Liew et al. (2005).Thus, in each of the three locations the model performance (C1 and C2) deteriorates in verification period.
The performance decrease is much more evident for the C1 model.This may be because the model computes melt as a function of only temperature and a melt parameter, whereas C2 also incorporates EI and V S .
In the application for Deadman Hill location, fig.( 14), the two models have similar performance in most years, but in several years (i.e. 2004-2006) the C1 model (classical temperature index) under-predicts SWE, potentially suggesting a stronger role of shortwave radiation in providing melt energy for those years.The applications in the Joe Wright and Hourglass sites, figures ( 15) and ( 16) respectively, show higher sensitivity in the C2 model, which generally has a higher peak SWE.However, for any of these applications, all parameters related to snow accumulation and melt have been optimized, so it is not possible to determine whether differences between the performance of the melt models is related primarily to the model structure or to the parameter combination selected in the optimization algorithm.

Test 2: SWE hourly simulations: models intercomparison
In this application the three snow melt components (C1, C2 and C3) were calibrated and verified against measured data.This can be performed in Joe Wright station, where air temperature, rainfall and snow water equivalent at hourly time step were available.The NewAge OMS component SWRB, Formetta et al. (2013b), was used to estimate incoming solar radiation time series input of C3 component.
In this case, the calibration period was 2008, and the verification period was from 2009 to 2013.
Table 16 presents the optimal parameter set for each component (C1, C2 and C3).For the T melt parameters, it identifies optimal negative values for model C1 and C2 and optimal positive value for model C3.The values are in line with the value founded in Kokkonen et al. (2006).The freezing coefficient, α f , assumes an optimal value between 0.0085-0.0099[mm C −1 h −1 ] for the three models.The coefficients for snow (α s ) and rainfall (α r ) precipitation are in line with the range estimated in literature (Forland et al. (1996) and Rubel and Hantel (1999)).

400
Table 17 present the indices of goodness of fit for calibration and verification periods of each model.Finally, figure 17 presents the comparison between simulated and observed snow water equivalent for the three models.The gray dots are the measured data, and green, red and blue solid lines are 405 the modeled data by using C1, C2 and C3 component respectively.
From the plot in fig.( 17) and the analysis of tab.( 17), it is clear that the three models are able to capture the variation in time of the snow water equivalent.Moreover, similar 410 to the daily time step application, the three models deteriorate their performance in the verification period.C2 is the best in preserving the goodness of fit, whereas C1 and C3 are the models that better capture the snow water equivalent for the calibration period and for 2012 portion of the time  17).The plot of measured snow water equivalent data shows how difficult it is to collect accurate measurements at an hourly 425 time-step: instability of the signal is evident in 2011 and 2013.Models, even simple as presented in this paper, could provide help in identifying data errors.

Test 3: Distributed application of SWE-C
The aim of this application is to show that SWE-C compo-430 nent is able to produce snow melt maps.Moreover the full integration in the NewAge JGrass system allows the user to immediately visualize and manage output raster maps.The SWE-C model is tested in distributed mode for the entire Cache la Poudre basin.The simulation period was between 1 435 October 2008 and 1 October 2009.Daily rainfall and temperature raster maps were computed using the detrended kriging algorithm.In this case three SNOTEL and three COOP meteorological stations were used, tab.( 11).The mean values of the three optimal parameters sets identified in the first The stable version of the model used in this paper will be available under GPL version 3 license at: http://code.google.com/p/jgrasstools/.The research version used in this paper is available on a GITHUB repository.
2103 SWE-C integration in the JGrass-NewAGE system of the data and setup of the simulations Three applications were performed: simulation of SWE at daily (Test n.1) and hourly (Test n.2) time step and model application in distributed mode (Test n.3).Calibration of SWE-C was conducted within the NewAge-JGrass system 320 shown in fig.( 415 series shown in fig.(17).In the last event (2013) in fig.(17),C1 and C3 models overestimate peak snow water equivalent, whereas C2 has a stronger performance.The timing of complete melting of the snow is well simulated by the three models, with Cazorzi and Dalla Fontana's model (C2) performing 420 the best of the three in the verification period Finally, hourly SWE data as available on the website http://www.wcc.nrcs.usda.gov/snow/are shown in fig.(

440
model test (section 4.3).Simulated SWE distributions for select dates are presented in fig.(18).Snow water equivalent maps were plotted for each month starting from 1 December 2008 to 1 April 2009, and six classes of snow water equivalent value are plotted for each month fig.(18).The model is able to simulate the spatial and temporal variability of the SWE during the hydrological year.The dynamics of the snow accumulation and melt are consistent with model structure and expected seasonal pattern of snow, where the snow accumulation increases with elevation and peaks in spring.Fu-450 ture work will compare these simulated SWE patterns to a more complex distributed model simulation of SWE and to satellite-derived snow cover patterns.5ConclusionsThis paper presents a parsimonious snow water equivalent model based on water and ice mass balance.The model simulates snow melt using one of three separate temperaturebased formulations where melt rates are a function of either temperature only or both temperature and solar radiation.The model is integrated into the NewAge-JGrass hydrological model as an OMS3 component, and for this reason it can make use of all the OMS3 components of the system: GIS based visualization, automatic calibration algorithm, and verification packages.All these components are applied and verified at three SNOTEL stations located in the Cache la Poudre river basin (Colorado, USA), and the model has good performance for both daily and hourly time steps, although model performance degrades from calibration to verification periods.Finally, the model is applied in distributed mode to simulate spatial patterns of SWE across the basin.Modeling snow water equivalent patterns in a distributed mode provides the possibility to compare them with more physically based snow models and the option to verify them with snow water equivalent remote sensing data.Future research will address problems related to modified temperature index snow water equivalent models such as transferability of parameter values to new locations and time periods, over-parameterization, comparison with physically based snow models, and verification of how well simulated snow cover spatial patterns reproduce spatial and temporal variability of the snowpack.

Fig. 11 .
Fig. 11.The NewAge System showing all the modeling components.Starting from the top: the uDig GIS, the meteorological data interpolation tools, energy balance, evapotranspiration, runoff production-routing and snow water equivalent.The user can select and connect different components and use automatic calibration algorithms (at the bottom) to optimize model parameters.

Fig. 12 .
Fig. 12.The SWE-C integration in the NewAge System showing connections with the short wave radiation component and kriging interpolation algorithm.Connection with the Particle Swarm Optimization algorithm is in red dashed line.

Fig. 13 .
Fig. 13.Cache la Poudre river basin digital elevation model.Elevations are in meters.

Fig. 14 .
Fig. 14.Calibration and verification model results for stationspecific calibration test in Deadman Hill station: the gray dots represent the measured SWE, the solid red line represents the model C1 (classical temperature index model) and the blue solid line represents the model C2 (Cazorzi and Dalla Fontana).

Fig. 15 .
Fig. 15.Calibration and verification model results for stationspecific calibration test in Joe Wright: the gray dots represent the measured SWE, the solid red line represents the model C1 (classical temperature index model) and the blue solid line represents the model C2 (Cazorzi and Dalla Fontana).

Fig. 16 .
Fig. 16.Calibration and verification model results for stationspecific calibration test in Hourglass station: the gray dots represent the measured SWE, the solid red line represents the model C1 (classical temperature index model) and the blue solid line represents the model C2 (Cazorzi and Dalla Fontana).

Fig. 17 .
Fig. 17.SWE-C application with hourly time time step in Joe Wright station.The gray dots are the measured data, and green, red and blue solid lines are the modeled data by using C1, C2 and C3 component respectively.
).Six meteorological stations have precipitation and temperature data in this river basin.These stations are presented in fig.13, and tab.11 shows their main features.Hourglass, Deadman Hill and Joe Wright are part of the Natural Resource Conservation Service Snow Telemetry (SNOTEL)

Table 11 .
Meteorological stations used in test simulations for the Cache la Poudre river basin.

Table 12 .
Optimal parameter values estimated for each of the three SNOTEL stations for the C1 model at a daily time step.

Table 13 .
Optimal parameter values estimated for each of the three SNOTEL stations for the C2 model at a daily time step.E is in W m −2 .

Table 16 .
Optimal parameter values estimated Joe Wright station for C1, C2 and C3 model.E is in W m −2 .

Table 17 .
Goodness of fit values for the calibration (top) and verification periods (bottom) for the site-specific calibrations at Joe Wright and for the models C1, C2 and C3 at a hourly time step.