The impact of different elevation steps on simulation of snow covered area and the resulting runoff variance

This study analyses the impact of vertical model discretisation on modelling snow covered area and the consequential effects on runoff formation of the semi-distributed water balance model HQsim. Therefore, the parameters relevant for snow modelling are varied within the frame of a uniformly distributed Monte Carlo Simulation (MCS). Since the model is based on the hydrological response unit (HRU) approach, the effect of building the HRUs with different elevation steps (250 m and 500 m) is tested for two alpine catchments. In total 5000 parameter combinations were generated for simulation. The results of modelled snow covered area were compared with thirty MODIS (Moderate Resolution Imaging Spectroradiometer) snow cover maps for the melting periods in 2003–2011. Based on a contingency table the comparisons were evaluated by different skill measures. Finally, the pareto optimal parameter settings of each skill measure were detected. Evaluation of runoff variability within the MCS and their pareto optimal runs show reduced variances of model output resulting from an improved simulation of the snow covered area.


Introduction
Snow represents an important component of the hydrological cycle.Especially for mountain regions the snow influences runoff generation in several ways.On one side, precipitation, stored as snow in winter, contributes essentially to runoff generation in spring and early summer.On the other side, snow has the potential to affect flood events.By storing melt water and rain, unsaturated snow acts as a buffer reducing the runoff (Schöber et al., 2012).In contrast, when rain and snowmelt add up, the runoff may increase rapidly (Weingartner et al., 2003;Woo, 2005).Hence, the estimation of snow is a considerable factor for rainfall-runoff modelling, especially for flood forecasting in alpine catchments.
Remotely sensed snow cover data provide a basis for calibrating conceptual hydrological models and improve the snow estimates of models without significant loss in runoff performance.Especially snow cover products of MODIS (Moderate Resolution Imaging Spectroradiometer) offer an appealing source because of their almost daily availability (Hall et al., 2002;Parajka and Blöschl, 2008a).Parajka and Blöschl (2008b) tested MODIS snow cover data in addition to runoff for calibrating a hydrological model.They found that the runoff performance of the calibration period was slightly poorer when calibrating on runoff and snow model efficiency than by runoff calibration only.In contrast, by using the results of combined calibration there was an improvement of runoff performance during the verification period.Hence, the calibration of the snow covered area is a vital tool for improving model performance, especially in alpine catchments.
This study was done in frame of the project: "HoPI" ("Hochwasserprognose für den Tiroler Inn") -flood forecasting system of the Tyrolean River Inn.The hybrid forecasting system consists of several models.Tributary catchments of the river Inn are modelled with the "snow and glacier melt model (SES)" for glacierized (Schöber et al., 2010) and HQsim for unglacierized areas (Achleitner et al., 2012).The runoff outputs of these tributary models are used as input for the hydraulic model for the River Inn.Especially for the flood forecasting of melting periods, a hydrological model, well calibrated regarding the snow aspects, could be able to improve the predictions on the longer term.Thereby an improvement of predictions is not merely a more realistic simulation, but rather a more certain forecast value and parameter setting, respectively.For overparameterized environmental models, as HQsim is, there exist several settings which are able to give good fits to data.This status is called equifinality (Beven, 2006) and among others it results in the uncertainty of the parameter values.Furthermore, in this context a more certain forecast value means that the runoff estimators of the acceptable parameter settings do not differ widely in one time step.For reduction of parameter uncertainty and prediction uncertainty, respectively, the knowledge of understanding of the model behaviour and separating of non-influential model parameters is prerequisite.Among other sensitivity analysis methods the principle of variance-based methods is to compare the total model output variance with the output variance by keeping one parameter fixed.Thereby the measure of parameter sensitivity is the degree of output variance reduction of each parameter (Saltelli et al., 2008).Based on the problem of equifinality and the principle of variance-based sensitivity methods, the certainty of the forecasted value can be described with the variance of all acceptable simulations.That means that the predicted simulation values with a minor variance are more certain than simulation values with a major variance.
The study analyses the applicability of MODIS-data for calibrating the parameters, relevant for snow modelling, of HQsim.Of particular interest are (i) the parameters distribution (ii) the variance of predicted runoff values by calibrating the snow model and (iii) the effect of vertical model discretisation thereon.This paper starts with a short presentation of the study area and the applied hydrological model in Sect. 2. Further the used methods for parameter variation and comparison of simulation versus remote sensing data are described therein.Section 3 contains the analysis of parameter uncertainty and Sect. 4 its impact on runoff simulation.

Study area and set-up of the hydrological model
The study is performed for the catchments Brandenberger Ache (280 km 2 , 508-2250 m a.s.l.) and Brixentaler Ache (330 km 2 , 500-2466 m a.s.l.) located in the north-eastern part of Tyrol (Fig. 1).Both tributary catchments of the river Inn were simulated with HQsim using a large range of parameter settings and hourly time series of meteorological inputs, rainfall and temperature, for the period 1996-2011.
The applied model HQsim (Kleindienst, 1996) is a semidistributed water balance model, based on the hydrological response unit (HRU) approach.For this study, we built the HRUs of each catchment based on elevation steps of 250 m (M250) and 500 m (M500).Each HRU represents similar characteristics with respect to exposition, soil, vegetation and elevation.The precipitation input of the two catchments was spatially interpolated on a 5 × 5 km grid.For each HRU the rain sum is the result of an intersection with these grid cells.The state of precipitation, whether it is falling as rain, rainsnow mixture or snow, depends on the air temperature, which is calculated for the mean altitude of each HRU.Measurements of the air temperature from different elevation steps are the input for the calculation of base temperature and according lapse rates (Rinderer et al., 2008).For determining the snowmelt a modified degree-day factor approach is applied.Therein, exposition and slope of the HRU and consequentially the angle of incidence of the direct shortwave radiation are included.Furthermore, the concept of "cold content" (Braun, 1985) defined as the sum of negative temperature accumulated over previous days is used for calculating the delayed onset of snowmelt, when air temperature rises above the melting point.The drainage of the snowpack is calculated by the Snow Compaction Relationship after Bertle (Bertle, 1966;Kleindienst, 1996;Knauf, 1980).More details of the hydrological model e.g. a list of the main parameters can be found in Achleitner et al. (2009Achleitner et al. ( , 2012)).
The hypsographic curve of the catchments (Fig. 1) shows that the largest part of the Brandenberger Ache is located between 1000 m and 1500 m.This has consequences for the area distribution of its models M250 and M500.About 70 % of the total area of the model M250 is located around the elevation steps of 1000 m-1250 m and 1250 m-1500 m.Even 65 % of total area of model M500 is located at the elevation step of 1250 m-1750 m.Furthermore, consequences of the vertical discretisation are the number of HRUs and their mean area.The models of Brandenberger Ache consist of 25 HRUs with a mean area of 12 km 2 (M250) and 14 HRUs with a mean area of 21 km 2 (M500), respectively.In contrast, the hypsographic curves of the catchment and the models of the Brixentaler Ache rise more steadily.The models consist of 54 HRUs with a mean area of 6 km 2 (M250) and 39 HRUs with a mean area of 8 km 2 .

Parameter variation using Monte Carlo Simulation
Considering physical characteristics, we defined feasible ranges of snow controlling parameters (Table 1).First, these are parameters to control the precipitation input to the model.Above or below the upper or lower temperature of snow/rain transition (tsrmin/tsrmax) the precipitation is falling totally as rain or snow.Between them, precipitation is calculated as a snow-rain-mixture.The snow correction factor (snowcorr) considers the systematic error of measurements causing most time underestimation of solid precipitation.Second, there are some parameters to describe the snow pack.In order to restrict the temperature of the snow pack there are two parameters: the lower limit of the snow temperature (sntmin) below which the temperature cannot decrease, and the maximum days (sntmem), which define how long the air temperature is stored for calculating the snow temperature of the next time step.The snow conditions are further described by the content of liquid water within the snow pack, which is limited by the maximal liquid content as an upper bound (maxliq).Third, the remaining parameters of the snow model regulate the snow melt.The constant "groundmelt" regards the snowmelt conditioned by heat influx of the underlying ground.The snowmelt induced on a degree-day basis can be adjusted with the minimum (snmfmin) and maximum (snmfmax) degree-day factor.The "albedo" describes the mean albedo of the surrounding area of the HRU.For the drainage of the snowpack calculated by the Snow Compaction Relationship after Bertle (1966), a critical density (critdens) is defined below which no liquid water can be stored in the snow pack.
For this study, only the snow controlling parameters mentioned above were varied separately by a uniformly distributed Monte Carlo Simulation (MCS).Parameter dependencies, such as the minimum degree-day factor must be lower than the maximum degree-day factor, were considered.Similar runs were avoided by applying the rule that every parameter must change more than ±5 % compared to any previous run.Therefore, the full parameter range was defined as one hundred per cent.A detailed description of setting the parameter values with percentage variation is given in Achleitner et al. (2009).In total, 5000 parameter combinations were generated for running HQsim in each catchment with both HRU discretisation (M250 and M500).

Snow cover estimation: simulation versus remote sensing
We used MODIS data with a spatial resolution of 500 m and a daily temporal resolution (Hall et al., 2002) for comparing our simulated snow covered area with remote sensing data.The MODIS instruments are installed on two satellites (TERRA and AQUA) which scan the same surface with a time shift of three hours.Global snow cover products, derived from MODIS observations, are available from the Distributed Active Archive Center of the National Snow and Ice Data Center (NSIDC, www.nsidc.org).For the snow melt periods of 2003-2011 daily snow cover products of both satellites were used.We selected one snow map with the least cloud cover interference for each month during all melting periods.Altogether 30 dates were selected.Figure 2 illustrates the procedure from the selected snow maps through the contingency tables to the statistical skill measures.Referring to Parajka and Blöschl (2008a) we combined the data of the Aqua and the Terra satellite for a first step cloud reduction.In a second step, the remaining cloudy cells were replaced with the value of the majority of eight adjacent cells.
To compare simulated results with MODIS data a spatial intersection of the MODIS grid and the HRUs was used.Therefore, the subareas of each HRU contributing to one MODISpixel were calculated.Per definition, a simulated pixel counts as snow covered, if the sum of the snow covered area is equal or higher than 50 %.For analysing the correspondence of modelled versus observed values the grids were evaluated by different skill measures based on a contingency table (Fig. 2).Five skill measures, recently presented in this context by Zappa (2008), were applied with n 11 = the sum of correctly predicted snow covered cells, n 00 = the sum of correctly predicted snow free cells, n 10 = the sum of false predicted snow covered cells, n 01 = the sum of false predicted snow free cells, n 1x = the sum of modelled snow covered cells, n 0x = the sum of modelled snow free cells, n x1 = the sum of observed snow covered cells, n x0 = the sum of observed snow free cells, n xx = the total sum of cells: the Accuracy (ACC), with a value of 1 for perfect correspondence of modelled versus measured values and 0 for no skill: The systematic error (BIAS) has its optimum at a value of 1.
Values below 1 indicate underestimation and above 1 overestimation: The False Alarm Ratio (FAR), the fraction of false simulated pixels has its optimum at 0: The Critical Success Index (CSI), the number of correct positive forecasts divided by the difference of all pixels and "hits" (CSI) ranges between 0 and 1 with its optimum at 1: And the Heidke Skill Score (HSS) with its score range between −1 and 1 and the best agreement at 1. (5) To obtain the best parameter settings, we first calculated every skill measure for each MODIS snow cover dataset for all simulation runs.Out of these, we identified the pareto optimal parameter settings (Gupta et al., 2005) of each of the five skill measures.That means that no parameter setting will produce a better skill measure in one event without a degradation of another event.Finally, we compared the five lists of pareto optimal parameter settings of the skill measures and selected those with the highest frequency of occurrence.The result of this procedure was a selection of parameter settings which are found in each list of every skill score.

Analysis of parameter uncertainty and its impact on runoff simulation
The prior uniform distribution of the parameters results in a posterior distribution after the pareto selection procedure.
For the models of Brandenberger Ache and Brixentaler Ache the posterior parameter distribution is illustrated in Fig. 3. Therein some of the distributions differ significantly from the initial one.The most prominent ones are the snow correction factor and the maximum value for the melt function, where the first one tends in all four cases to a very right skewed distribution.On the other hand, the melt function maximum has got a very left skewed distribution.Its counterpart, the melt function minimum reaches a symmetrical distribution around value 4 for Brixentaler Ache.Still these two parameters are special cases, due to their link to each other in the generation of parameters.Small values for the melt function maximum were eliminated due to the constraint that the minimum degree-day factor must be lower than the maximum degreeday factor.All the other parameters do not have such a clear distribution, but tendencies are detectable.For example, the maximum snow-rain transition tends to values between 0 • C and 2 • C or the days of temperature memory has got an accumulation between 0 and 10 days.
The effect on variances of simulated runoff for the melting period of 2010 is given in Fig. 4. The Fig. 4a shows the gauge "Mariathal" (Brandenberger Ache) and Fig. 4b illustrates the gauge "Bruckhäusl" (Brixentaler Ache), respectively.During the melting period, between the 12 April and the 25 May, the runoff regime is seasonal adjusted accordingly.At the end of May 2010 a cold front reached the northern part of Tyrol.Above 1000 m a.s.l.cool temperatures led to solid precipitation and a thin snow cover.Due to large precipitation amounts floods occurred at Brandenberger Ache (> 10 yr return period) and Brixentaler Ache (> 5 yr return period) (Gattermayr, 2010a, b).The Fig. 4c and d illustrate the corresponding normalized variances of the simulated runoff of the models M500 and the plots Fig. 4e and f of the models M250, respectively.There are obvious differences between the models M500 and M250 of Brandenberger Ache.Especially in the period of the 12 April and the 25 May, the variances are reduced significantly.This leads to the assumption that the model M250 did not accumulate as much snow as the model M500.The plots of normalized variances of the Brixentaler Ache look very similar.Most visible differences between the discretisations of 250 m and 500 m are detectable in the period of the 2 May and 11 June.In this period the variances of the model M250 is slightly lower than of the model M500.In all four cases the graphs shows a significant reduction of runoff variance due to acceptable snow simulations.The most notable reduction is seen at the models M500 of the Brandenberger Ache in the melting period between the 12 April and the 22 May.After this date the normalized variances of the pareto optimal parameter settings are so low that they are not shown in the graph.For the model M250 the reduction of variance is not as obvious as for the model M500.Especially in the beginning of the melt season, between the 12 April and the 2 May, the runoff variances of the acceptable snow simulations are only slightly lower than of all MCS.For the graphs of the Brixentaler Ache the situation is related but not as pronounced as Brandenberger Ache.In both catchments the runoff variances of the models M500 are more significantly reduced by the pareto optimal parameter settings than the variances of the models M250.This effect can be seen in the runoff graphs of the event in June 2010 in Fig. 4g-j.Therefore we do not take into account the relation to the discharge measurements consciously, because we do not vary the full parameter set effecting runoff simulations in this study.Beyond that the reduction of runoff variability of the pareto optimal settings are of main interest so the range of simulated runoff of the whole MCS is confronted with those of their pareto optimal simulations.The graphs of the Brandenberger Ache (Fig. 4g and i) show clearly, that the range of the models M500 is much wider and tend to simulate to higher runoff than those of the models M250.In contrast, the range of the pareto optimal runs of the models M500 is closer than the ranges of the models M250, especially afterwards the 2 June.For the models of the Brixentaler Ache (see Fig. 4h and j) the effects are not as obvious as for the models of Brandenberger Ache.Most significant differences between the models M250 and M500 of the Brixentaler Ache occur between the 1 June to the 3 June and between the 4 June and 8 June.In the first period the range of the pareto optimal runs of the models M250 is slightly closer than of the models M500.In contrast, in the second period the ranges of the pareto optimal runs of the models M500 exhibit a less variability in the diurnal cycle than the models M250.In this period the effect of simulate the snow covered area acceptable shows a significant influence on the runoff of the models M500.

Benefits and limitations using remote sensing data for improved model calibration
From the MCS the pareto selection procedure led to varying frequencies of parameter settings of acceptable simulations.
For Brandenberger Ache we selected 121 parameter settings  for the model M250 and only 19 for the model M500.The Brixentaler Ache simulated the snow covered area with much more acceptable parameter settings.1225 combinations for the model M250 and still 730 settings for M500 were found.There are 84 parameter combinations for model M250 and 12 settings for M500 with acceptable simulations in both catchments.This shows that (i) applying a vertical discretisation of 250 m results in more parameter settings found to be pareto optimal than using the 500 m discretisation and (ii) the amount of pareto optimal parameter settings for Brixentaler Ache is at least more than ten times higher than for Brandenberger Ache.This can be explained by: 1.In general the HRUs of M250 cover a smaller area than HRUs of M500.Hence, because of the intersection between the MODIS grid and the HRUs, a HRU of M250 takes part in less simulated MODIS cells than a HRU of M500.This results in a less relative impact of false modelled HRUs of M250 within the simulated MODIS cells.
2. About 60 % of the area of Brandenberger Ache is located between 1000 m a.s.l. and 1500 m a.s.l.In particular in melting periods this elevation range is a sensitive zone for the extent of snow cover.Especially for modelling snow covered area with the model M500 and its highly extended HRUs there are just a few parameter settings which are able to represent the snow pattern correctly.In contrast, for Brixentaler Ache, with an approximately linear increasing hypsographic curve, the pattern of snow covered area can be modelled with more parameter settings which results in a higher variability of runoff simulations.
In a holistic view, MODIS-data appertain for calibrating the snow model of the hydrological model HQsim.The selection of the parameter settings by acceptable simulations has shown that the comparison with snow maps reduces the count of settings significantly.Moreover, the variance of the simulated runoff decreases in comparison with all results of the MCS.

Conclusions
The aim of this study is to analyse the impact of the vertical model discretisation on the variability of runoff, which are introduced by the correctness of the modelled snow covered area.Thereby the focus is set on the posterior parameter distribution, the variance of the predicted runoff values and finally the effect of the vertical discretisation thereon.Therefore, a uniform distributed MCS of the snow relevant parameters was computed and their snow covered area results were compared with 30 snow maps altogether.The results show that the number of parameter settings for well snow covered area simulation can be reduced by calibration with MODIS-Data.In particular for catchments with no steady increase of area with elevation the parameter settings will be reduced significantly.Thereby the vertical discretisation has a major influence because of covering critical elevation steps with at worst only one elevation zone (Fig. 1).For catchments with an approximately linear increasing hypsographic curve the vertical discretisation of the model seems to have not such an influence like at Brandenberger Ache.Also, at the catchment of Brixentaler Ache the count of parameter settings between the model M250 and M500 do not differ as much as at Brandenberger Ache.The results of the calculated output variance show similar tendencies.The most significant reduction of variance is detectable at the model M500 of Brandenberger Ache followed by the model M500 of Brixentaler Ache.Hence, MODIS snow maps are a very helpful tool for calibrating the snow relevant parameters of a hydrological model.Especially the variance of the model output in snow melting periods can be reduced considerably and the predicted values are more certain in this way.Finally the variance of the simulations can be explained by the different snow depths which lead to variability in the stored snow water equivalent of the system.

Fig. 1 .
Fig. 1.Overview on the location of the modelled catchments and their hypsographic curves.

Fig. 2 .
Fig. 2. Procedure of the comparison of simulations versus MODIS snow maps (map colours: grey: cloud, white: snow, black: no snow).

Fig. 3 .
Fig. 3. Histograms of the distribution of the Pareto optimal parameter settings.A description of the parameters can be found in Sect.2.2.

Fig. 4 .
Fig. 4. Hydropgraphs of the melting period 2010 for (a) Brandenberger Ache (gauge: "Mariathal") and (b) Brixentaler Ache (gauge:"Bruckhäusl").The second row illustrates the normalized variances of simulated model output for the models M500 and the third row for the models M250, respectively.The two rows below illustrates the simulated runoff of the event in June 2010 (grey: entire Monte Carlo Simulation, red: optimal pareto settings).

Table 1 .
Varied snow relevant parameter of HQsim by Monte Carlo Simulation.A description of parameters can be found in Sect.2.2.