The Cache la Poudr river basin snow water equivalent modeling with NewAge-JGrass

Introduction Conclusions References


Introduction
The physically based distributed approach is the best way to simulate the snowpack evolution.This solution has reached maturity and was pursued successfully with many recents models including CROCUS (Brun et al., 1992), Alpine3D (Lehning et al., 2006), GEOtop (Rigon et al., 2006;Zanotti et al., 2004;Endrizzi, 2007;Dall'Amico et al., 2011), ISNOBAL (Marks et al., 1999), UEB (Tarboton and Luce, 1996).These models often implement, besides the core energy budget, ancillary modeling of blowing snow, and other features that are required to reproduce the full set of thermodynamic quantities representative of snowpack state.However, performing the snow budget and modeling its complete variability is not always necessary and requested.In many situations, for instance where the prognostic significative quantity is just the total snow water equivalent in a sub-catchment, more simple models can work better.Besides, realtime modeling with data assimilation and parameter calibration requires that a whole forecasting cycle is obtained in few minutes for an entire day, in order to proceed with all the appropriate operations and testing.In any case, a best practice is to compare the most complete models with the simplest ones in order to assess the degree of complexity that is required for any task.The ancestor of all these simple models, is the SRM model by Martinec (1975) which was implemented several times and applied 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 an average of the daily temperature and an estimate of the catchment area covered by snow which is called snow water depletion curve.These are tricky to determine, but possible to be detected by satellites.Therefore the model was largely used together with remote sensing data.In this paper we implement one of the minimalist SWE models based on the idea, partially investigated in Zanotti et al. (2004).Once a good estimate of radiation is available, good spatially distributed estimates of the snow-water equivalent can be obtained.Introduction

Conclusions References
Tables Figures

Back Close
Full  et al. (1994) andBrubaker et al. (1996) introduced simple SWE modelling based on the use of the radiation budget.However, in this dissertation we use the formulation of the problem developed by Cazorzi and Dalla Fontana (1996), since it was based simply on the estimate of the direct solar radiation, rather than the total net radiation, which is more rarely used and more difficult to obtain.

Kustas
The SRM parameters, in the original philosophy adopted by its authors, were not calibrated or optimised by historical data.They could be either derived from measurements or estimated by hydrological judgement taking into account the basin characteristics, physical laws and theoretical or empirical relations.In many studies, this hypothesis was relaxed, and we adopt a completely opposite strategy, in which we use all the available data to assess the model's parameters.Therefore, we make use of data measured at stations and use the particle-swarm optimizer, proposed by Kennedy and Eberhart (1995), to obtain the parameters of the model which can be, eventually, studied for detecting regularities, and gaining insights about the phenomena studied.
Another novelty of our model, NewAge-SWE, is that it does not come alone, but as part of a larger system, called NewAge-JGrass (Formetta et al., 2012 andFormetta et al., 2013, which includes several modeling components, briefly described in the next section of the paper.NewAge-JGrass, in turn, is based on the Object Modeling System version 3 (OMS3) modeling infrastructure (David et al., 2002).The OMS3 infrastructure promotes the modern concepts of object-oriented programming applied to hydrological modeling.Using it, any part of a model can be deployed as a "component" which can eventually be connected, just before the run-time, to other components with a scripting language to provide a modeling composition suitable to solve the problem under examination, and to be compared to other modeling compositions obtained with alternative components.Introduction

Conclusions References
Tables Figures

Back Close
Full The component presented in this paper are part of the NewAge-JGrass system (Formetta et al., 2011).It includes components dealing with the various hydrological processes such as estimation of: the space-time structure of precipitation, (krigings) shortwave and long wave radiation balance, (SWRB and LWRB) evapotranspiration, (Priestley-Taylor or Penman-Monteith) runoff production, (Hymod) aggregation and propagation of flows in channel (Routing) Finally, it also includes three different automatic calibration algorithms: -Particle Swarm Optimization component (Eberhart and Shi, 2001) -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 evaluation is the hillslope.Each hillslope drains into a single associated link rather than cells or pixels.The model's physics requires the 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 implements algorithms that take into account shadows and complex topography.Shortwave radiation under generic sky conditions (all-sky) is Introduction

Conclusions References
Tables Figures
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 component's behavior is eventually tested with the use of the component NewAge-V (Verification) which provides some of the classical indices of goodness of fit such as: Nash-Sutcliffe, Percentage bias, Index of agreement and Kling Gupta efficiency.
Finally, every component can be connected, parameterized, and executed either using the OMS3 console (OMS 3.1) or the OMS3 scripting mode within the uDig Spatial Toolbox (http://code.google.com/p/jgrasstools/).A custom hydrological model obtained by joining several components is usually called "modeling solution".It can be instantiated, initialized and connected in a sequence.In this way the modeler can build a custom hydrological model and solution by selecting among the components to simulate those useful to solve the hydrological problem under analysis.The model composition obtained, once executed, will use the OMS3 implicit parallelism to improve the computational efficiency in multicore or multiprocessor machines.

The NewAge-SWE component's equations
The snow melting model is based on a modificated approach presented in Kokkonen et al. (2006).The main novelties are: melt formulation: in the presented version melt is a function not only of temperature but also of a radiative term as presented in Cazorzi and Dalla Fontana (1996); in Kavetski et al. (2006), it could generate extremely non-smooth parameter surface during the automatic calibration procedure.In order to avoid this problem a smoother for thresholds is used; model structure: as integrating part of the OMS3 NewAge model infrastructure, the presented snow melt component can be applied at the hillslope scale and can utilize all the tools such as meteorological data interpolation, automatic calibration and radiation balance algorithms.
In the next subsection the main algorithms of the model are described.

Mass balance
The snowpack mass balance is written as follows.For the water equivalent of ice it results: and for liquid water (M w [L]) in the snowpack it results: Equation ( 1) represents the variation in the time of the ice in the snowpack which is equal to the algebraic sum of the snowfall, P s , freezing, F , and melting, M (all expressed as snow water equivalent).Accordingly Eq. ( 2) the variation in time of the liquid water in the snowpack is equal to the algebraic sum of the rainfall, P r , freezing, F , and melting, M. If liquid water M w exceeds liquid water-retention capacity of the snowpack (M max [mm]), the surplus becomes snowmelt discharge q m [L T −1 ].The liquid water retention capacity of a snowpack is related to the ice content by a linear relationship, Eq. ( 3)

Conclusions References
Tables Figures

Back Close
Full Differently from Kokkonen et al. (2006), the time step to be used in these two coupled equations is not necessarily daily: but this is made dependent on the interval at which the fluxes in the second member of Eqs. ( 1) and ( 2) are made available.

The type of precipitation
The first hydrological process to be simulated is the discrimination between rainfall and snowfall considering that the two forms of precipitation appears as distinct in Eqs. ( 2) and ( 1).Usually only rain gauge measurement and air temperature are available.
A common procedure is to consider a threshold for the air temperature 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 As proposed in Kavetski et al. (2006) to avoid problems for parameter calibration, a smoother filter for thresholds is applied and the algorithm to discriminate between rainfall and snowfall can be described as follows: where: ] 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 region were presented (Forland et al., 1996;Rubel and Hantel, 1999;Michelson, 2004), in the model the two coefficients are considered parameters and therefore calibrated.Introduction

Conclusions References
Tables Figures

Back Close
Full

Snow melt fluxes
Based on the approach presented in Cazorzi and Dalla Fontana (1996) the melting process, Eq. ( 5), is a function of both shortwave radiation and air temperature.The two main differences in the presented model compared to Cazorzi and Dalla Fontana (1996) are: a new algorithm is used to compute the shortwave radiation (direct plus diffuse component) proposed by Corripio (2003) and integrated into NewAge-JGrass model Formetta et al. (2012) which accounts for the complex topography, shadows and the sky view factor Corripio (2002), and the cloud cover.The equation for the melt process is: where: is the energy index and V s [-] is the sky view factor.The energetic index is the potential energy accumulated over a given period at a certain point.To compute the energy index the shortwave energetic balance component implemented in NewAge Formetta et al. (2012) is used.The shortwave beam and diffuse solar radiation is accumulated for each pixel and the result is divided by the given period of the time.As presented in Cazorzi and Dalla Fontana (1996) five energetic index maps are computed starting from 21 December (winter solstice) to the middle of each month from February to June.During the night the snow melt is a function of the energetic index minimum value of the considered map, as presented in Cazorzi and Dalla Fontana (1996).

Freezing
The rate of freezing F that is compared in the mass budgets is linear related to the air temperature when the air temperature is less then the melting temperature, as

Conclusions References
Tables Figures

Back Close
Full presented in Eq. ( 6) where F [L T −1 ] is the freezing rate and α f [L • C −1 T −1 ] is the freezing degreeday(hourly) factor.If the model is used with daily time steps temperature is the mean daily temperature.If it is used at hourly scale, temperature is the mean hourly temperature.Accordingly the value of the parameter α f change values.

SWE-C integration in NewAge System
The SWE-C is perfectly integrated in the NewAge System as presented in Fig. 1.Firstly, it uses the meteorological interpolation algorithms: Krigings tools, for temperature and precipitation interpolation, and JAMI for the temperature interpolation.Like the interpolation algorithms, SWE-C is able to work at a raster and a point scale.Secondly it uses the NewAge short wave radiation component in order to estimate the maps of cumulated energy in different periods of the year as explained in the model equations section.This components is able to take into account complex topography, shadow, and clouds cover.Thirdly, the SWE-C outputs could be: raster maps or time-series (one for each hillslopes centroids) of snow water equivalent and snow melt.Those could be used by the rainfall-runoff components in order to model a river basin where the snow contribution is not negligible.Finally, the SWE-C component could be connected to the NewAge and OMS3 calibration algorithm in order to estimate the best model parameters values.Introduction

Conclusions References
Tables Figures

Back Close
Full

NewAge-SWE verification
The model is applied in the Cache la Poudre River basin.Three applications are presented in this section.Firstly, the model was applied point mode for three stations where snow water equivalent measurement were available the model was calibrated and verified.Secondly, simulations were performed in order to investigate how representative the optimal parameter sets are for each stations.Finally, the model is applied in the fully distributed mode: raster maps of the snow water equivalent over the entire basin are simulated.

Sites description
Model tests sites are located within the upper Cache la Poudre basin in the Rocky Mountains of northern Colorado and southern Wyoming, USA.This 2700 km 2 basin has elevations ranging from 1590-4125 m, with mean annual precipitation ranging from 330 mm at lower elevations to 1350 mm at the highest elevations.Six are the meteorological stations available on the river basin.They are presented in Fig. 2 and Table 1 shows their main features.Hourglass, Deadman Hill and Joe Wright belong to the Natural Resource Conservation Survey Snow Telemetry (SNOOTEL) meteorological stations.They provide data (precipitation, air temperature and SWE) at daily time step.For Hourglass station the data available start on 1 October 2008 and ends on 1 May 2012 (the first year is used as calibration period and the last 3 yr are used as validation period); for Joe Wright and Deadman Hill stations they go from 1 October 1999 to 1 October 2009 (the first year is used as calibration period and the last 9 yr are used as validation period).
Buckhorn were used for air temperature and precipitation interpolations in the fully distributed application of the snow melting and snow water equivalent component.

Test 1: model calibration and verification
As mentioned in the basin description there are three snow telemetering (SNOTEL) stations, Fig. 2: Hourglass, Joe Wright and Deadman Hill.Table 1 shows their main features.They provide daily rainfall, temperature, and snow water equivalent data.
For Hourglass station the available data starts on 1 October 2008 and ends on 1 May 2012 (the first year is used as calibration period and the last 3 yr are used as validation period); for the Joe Wright and Deadman Hill stations data goes from 1 October 1999 to 1 October 2009 (the first year is used as calibration period and the last 9 yr are used as validation period).
To calibrate the SWE-C the configuration of the NewAge-JGrass components shown in Fig. 1 was used.For this task, the the calibration algorithm Particle Swarm Optimization was used Kennedy and Eberhart (1995) and Eberhart and Shi (2001).
As objective function the Kling-Gupta Efficiency (KGE) presented in Gupta et al. (2009) was selected.
The model was verified for the three stations in two different ways.In a first approach a different optimal parameters set was estimated at each station and was used to simulate the validation period.The second method estimated the optimal parameters set in one station to model the simulation period in the other 2 stations and the procedure was repeated for each stations.For the Deadman Hill and Joe Wright stations the calibration period was the year 1999 and for Hourglass was the year 2008.
Theree classical GOF index are computed: Nash-Sutcliffe (NSE), Percentual Bias (PBIAS) and Index of Agreement (IOA).NSE values greater than 0.75 mean that the model can be considered "good" values between 0.75 and 0.36 are associated with a "satisfactory" model and values below 0.36 indicate "not a satisfactory" model.Looking at the hydrological mode classification, as presented in Stehr et al. (2008) and Van Liew et al. (2005), a model which presents an absolute PBIAS value less then 20 is

Conclusions References
Tables Figures

Back Close
Full considered "good", if the values are between 20 and 40 it is considered "satisfactory", and if it is greater than 40 the model is considered "not satisfactory".Table 3 shows, for the calibration period, at the top, and for entire simulation period, at the bottom, the indexes of goodness of fit for the three SNOTEL stations.
The model calibrated at each station and verificated by using the optimized parameter can be considered "good" in both calibration and validation periods even if the model performance in the validation period is slightly lesser.

Test 2: parameters representativeness
In order to investigate how representative a parameter set really is a number of simulations are performed.For the entire simulation period the optimal parameter set for the Deadman Hill station was used for estimating the other two stations and the GOF indexes were computed.The same methodology was also applied for the Hourglass and Joe Wright stations, respectively.The simulations results are presented in Table 4: the column "Optimal parameter set" specifies the station's parameter set used in the simulation.
As presented in Table 4, the model results are sensible to parameters variations.Even if the model for all the simulations performed can be classified as at least "satisfactory" for the NSE and PBIAS GOF's, this application emphasizes that the modeller has to pay attention to the parameters representativeness expecially at different locations.This becomes more inportant when the parameters are stictly related to measurement site features.For example, T m could depend on the elevation, aspect of the measurement site, α s and α r could be function of the measurement instrument, α m and α f could be connected to the causes related to the amount of energy collected at the site (sky view factor, vegetation, or antrophic occlusions)

Conclusions References
Tables Figures

Back Close
Full The SWE-C model is tested in distributed mode for the Poudre river.The simulation period was between 1 October 2008 and 1 October 2009.Daily rainfall and temperature raster maps were computed by using the detrended kriging algorithm.In this case three SNOTEL and three COOP meteorological stations were used.Table 1 shows their main features.
The mean values of the three optimal parameters set as presented in the previous section were used in this simulation.The results are presented in Fig. 6.Snow water equivalent maps were plotted for each month starting from 1 December 2008 to 1 April 2009.Six classes of snow water equivalent value are plotted for each month (Fig. 6).

Conclusions
In the paper a parsimonious snow melting and snow water equivalent model based on water and ice balance is presented.Here, the snow melt takes not only into account the temperature but also the energy received at the simulated point.The model is integrated into the NewAge-JGrass hydrological model as 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 validation packages.All these components are applied and verified at three SNOTEL stations located in the Cache la Poudre river basin (Colorado, US) providing satisfactory results at all sites.A second model application focuses on the parameter representativeness.It shows that extending optimal parameter set at some location decreases model performances especially when the parameters are strictly related to the climate and geomorphological features of the site.Finally, the distributed application in the Poudre river basin is presented.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.These tasks are left to further papers.Introduction

Conclusions References
Tables Figures

Back Close
Full  Full Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | form of precipitation determination: in the presented version the model does not work with a threshold separation based on air temperature because, as presented Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Mountain, Rustic and Virginia Dale belong to Service Cooperative Observer Program (COOP) meteorological stations.They only provide precipitation and air temperature.For the three stations the data available start on 1 October 2008 to 1 October 2009.Those data integrated the SNOOTEL stations measurements.They Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Formetta, G., Mantilla, R., Franceschi, S., Antonello, A., and Rigon, R.: The JGrass-NewAge system for forecasting and managing the hydrological budgets at the basin scale: models of flow generation and propagation/routing, Geosci.Model Dev., 4, 943-955, doi:10.5194/gmd-4-943-2011,2011.4450, 4451 Formetta, G., Rigon, R., Chávez, J. L., and David, O.: Modeling short wave solar Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

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

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

Fig. 4 .
Fig. 4. Validation model results in Joe Wright station: the gray dots represent the measured SWE and the solid black line represent the modelled SWE.

Fig. 4 .
Fig. 4. Validation model results in Joe Wright station: the gray dots represent the measured SWE and the solid black line represent the modelled SWE.

Fig. 5 .
Fig. 5. Validation model results in Hourglass station: the gray dots represent the measured SWE and the solid black line represent the modelled SWE.

Fig. 5 .
Fig. 5. Validation model results in Hourglass station: the gray dots represent the measured SWE and the solid black line represent the modelled SWE.

Table 1 .
List of the meteorological stations used in the simulations performed in the Cache la Poudre river basin.