the Creative Commons Attribution 4.0 License.

the Creative Commons Attribution 4.0 License.

# TOPMELT 1.0: a topography-based distribution function approach to snowmelt simulation for hydrological modelling at basin scale

### Mattia Zaramella

### Marco Borga

### Davide Zoccatelli

### Luca Carturan

Enhanced
temperature-index distributed models for snowpack simulation, incorporating
air temperature and a term for clear sky potential solar radiation, are
increasingly used to simulate the spatial variability of the snow water
equivalent. This paper presents a new snowpack model (termed TOPMELT) which
integrates an enhanced temperature-index model into the ICHYMOD
semi-distributed basin-scale hydrological model by exploiting a statistical
representation of the distribution of clear sky potential solar radiation.
This is obtained by discretizing the full spatial distribution of clear sky
potential solar radiation into a number of radiation classes. The computation
required to generate a spatially distributed water equivalent reduces to a
single calculation for each radiation class. This turns into a potentially
significant advantage when parameter sensitivity and uncertainty estimation
procedures are carried out. The radiation index may be also averaged in time
over given time periods. Thus, the model resembles a classical
temperature-index model when only one radiation class for each elevation band
and a temporal aggregation of 1 year is used, whereas it approximates a fully
distributed model by increasing the number of the radiation classes and
decreasing the temporal aggregation. TOPMELT is integrated within the
semi-distributed ICHYMOD model and is applied at an hourly time step over the
Aurino Basin (also known as the Ahr River) at San Giorgio (San Giorgio Aurino), a
614 km^{2} catchment in the Upper Adige River basin (eastern Alps, Italy) to
examine the sensitivity of the snowpack and runoff model results to the
spatial and temporal aggregation of the radiation fluxes. It is shown that
the spatial simulation of the snow water equivalent is strongly affected by
the aggregation scales. However, limited degradation of the snow simulations
is achieved when using 10 radiation classes and 4 weeks as spatial and
temporal aggregation scales respectively. Results highlight that the effects
of space–time aggregation of the solar radiation patterns on the runoff
response are scale dependent. They are minimal at the scale of the whole
Aurino Basin, while considerable impact is seen at a basin scale of
5 km^{2}.

Seasonal snow cover is important as storage and source of meltwater for human use, irrigation and hydropower production in many regions of the world. On the other hand, snow cover and meltwater can be a cause of disastrous natural hazards, such as floods and avalanches. Additionally, snow cover is a key factor in the weather and climate system, both regionally and globally (Armstrong et al. , 2008). Owing to society's strong need for updated information on snow conditions, snow accumulation and melt, models have been developed with a wide range of features. Approaches for snowpack computation range from empirical models (e.g. simple temperature-index models) to more-sophisticated physically based energy-balance models (Zappa et al., 2003; Vionnet et al., 2012; Essery et al., 2013; Magnusson et al., 2015; Avanzi et al. , 2016). Temperature-index models require only temperature as input and are based on the assumption of a linear relationship between this variable and melt rates, whereas energy-balance models are based on the computation of all relevant energy fluxes at the snowpack surface, and thus require extrapolation of numerous meteorological and surface input variables at the local scale (Jóhannesson et al., 1995; Cazorzi and Dalla Fontana, 1996; Hock, 1999; Hock and Holmgren, 2005; Anslow et al. , 2008; Formetta et al., 2014). Advances over the simple dependence of melt on air temperature by addition of radiation terms have been suggested in the last decades (Hock, 1999; Pellicciotti et al., 2005; Carturan et al., 2012), and are termed enhanced temperature-index (ETI) models here. In contrast to simple temperature-index models, where melt varies in space only as a function of elevation (given by temperature lapse rates), ETI models include a term for clear sky potential solar radiation. This term accounts for topographic effects (e.g. aspect, slope and shading) on the spatial distribution of melt, without the need for additional meteorological variables (e.g. global radiation and cloud data). ETI models have been found to provide a better representation of the spatial and temporal variability of melt controlled by solar radiation, when compared with simple temperature-index models, as reported by a number of authors (Cazorzi and Dalla Fontana, 1996; Hock, 1999; Pellicciotti et al., 2005; Carenzo et al., 2009, among others). Some of these approaches also better cope with the physical character of the melt process and provide a promising approach to modelling the snowpack at the catchment scale with fewer input data than energy-balance models, but allow for better model parameter transferability than standard temperature-index models (Carenzo et al., 2009).

In spite of their improved accuracy compared to simpler approaches, ETI models have been so far not integrated within lumped or semi-distributed, basin-scale hydrological modelling schemes, which are still frequently used to model sparsely gauged mountainous catchments. Integration of ETI snowpack models into lumped or semi-distributed hydrological models may have the potential to increase spatial transferability of calibrated snowpack model parameters for hydrological applications over ungauged mountainous basins, as shown by Comola et al. (2015). Another important implication of the stronger physical basis of the ETI model with respect to simpler “degree day” models is that it might be more appropriate for the study of the climate-change impact on melt regimes, as shown by Pellicciotti et al. (2005). Finally, increasing the accuracy of the modelled snow water equivalent may improve the outcomes of data assimilation procedures of remotely sensed snow cover information. In a few cases, semi-distributed models incorporating a spatial discretization in classes of elevation and aspect have been applied to represent the effect of exposition on snowmelt and adjust snowpack parameters accordingly (Klok et al., 2001; Konz and Seibert, 2010; Abudu et al., 2016). In other cases, a mean value of radiation has been used over the basin area in the same elevation band, modifying accordingly the melt parameters (Li and Williams, 2008). However, these types of tessellations have no allowance for representing the actual variations of radiation distribution over space and time, which is an important feature in ETI models (Pellicciotti et al., 2005).

This work describes a novel snowpack model (termed TOPMELT herewith), which integrates the ETI snowpack method originally developed in a spatially distributed way by Cazorzi and Dalla Fontana (1996) within a semi-distributed basin-scale hydrological model. In the model developed by Cazorzi and Dalla Fontana (1996), local snowmelt is computed by using a combined melt factor which is multiplied by a radiation index and positive air temperature. With TOPMELT, pixels with a similar radiation index and air temperature are identified by subdividing basin elevation bands into a number of radiation index classes. Then, the snowpack modelling is carried out for each class of radiation index and for each elevation band. This ensures achieving significant computational efficiency, which characterizes the temperature-index models, allowing at the same time for a stronger physical basis for ETI models. This is a potentially significant advantage when several model simulation runs should be carried out, such as in Monte Carlo-based parameter sensitivity and uncertainty estimation procedures.

The model accounts for the temporal variability of the radiation index by using local mean values of the index computed over given temporal aggregation intervals, ranging from 1 to several weeks. This means that a time-averaged solar radiation distribution is used over a given temporal interval, before substituting it with a new averaged distribution. With decreasing the updating interval, the accuracy of the model is expected to increase at the expense of the computational efficiency.

As the spatial distribution of clear sky solar radiation changes with time, a radiation class computed over two different periods may sample two different portions of the elevation band. This means that a pixel belonging to a certain class at a given time will belong to a different class at another time. TOPMELT incorporates a time-integration routine, which accounts for the temporal variability of the radiation index distribution, ensuring a consistent temporal simulation of the snowpack. Thus, TOPMELT permits full implementation of the ETI snowpack method, taking into account the seasonal evolution of the spatial distribution of solar radiation. Moreover, it provides a spatially continuous mapping of simulated snow water equivalent, in spite of the computationally efficient semi-distributed representation of basin-scale snowpack modelling. Depending on the number of radiation classes which are used in the model, the snowpack model makes use of solar radiation values which are spatially averaged over different areas. Effectively, the model resembles a classical temperature-index model when only one radiation class for each elevation band is used, whereas it approximates a fully distributed model with increasing the number of the radiation classes (and correspondingly decreasing the area corresponding to each class).

This paper describes in detail the structure of TOPMELT and of the
time-integration routine. The integration of TOPMELT within the ICHYMOD
hydrological model is also illustrated. Finally, results are reported from
the application of TOPMELT over the 614 km^{2} Aurino Basin at San Giorgio
in the Upper Adige River system (eastern Italian Alps). The case study is
exploited to (i) examine the sensitivity of the snowpack and runoff model
results to the temporal and spatial aggregation of the radiation fluxes, and
(ii) to identify suitable spatial and temporal aggregation intervals for
model simulation. The sensitivity analysis is performed on modelled snowpack
in terms of snow water equivalent and on the ensuing simulated runoff,
comparing the output from simulations performed at different aggregation
intervals with a reference represented by the finest aggregation levels.

In TOPMELT, the basin area is subdivided into elevation bands to account for
air temperature variability with elevation. Then, each elevation band is
subdivided into a number of radiation classes. This is carried out by
dividing each elevation band into a number *n*_{c} of equally
distributed radiation classes, where the *i*th class contains the band
sub-area corresponding to the *i*th percentile of the incident radiation
energy. Therefore, the model spatial domain is represented by *n*_{b}
elevation bands and by *n*_{c} radiation classes for each elevation
band. TOPMELT deals with separate snow and glacier melt: to account for the
presence of a glacier area associated to an energy class, each one of the
*n*_{b}×*n*_{c} model cells is characterized by the
corresponding fraction of glacier area. The spatial subdivisions controls the
balance between computational efficiency and model accuracy in the snowpack
model.

The following sections describe the main input and modules of the model, where an hourly temporal interval is used for model computations.

## 2.1 Clear sky potential radiation computation and derivation of radiation distributions

For the application of TOPMELT presented in this work, clear sky shortwave
solar radiation (W m^{−2}) is computed at each element of the digital
terrain model (DTM) by taking into account shadow and complex topography,
calculating the apparent sun motion (Swift , 1976; Lee, 1978; Oke, 1992) and the
intersection of radiation with topography (Dubayah et al., 1990; Ranzi and Rosso, 1991). Diffuse
radiation is computed by accounting for self-shading (by slope and aspect)
and occlusions produced by the visible horizon. Since the model uses
radiation values averaged over a given time interval, maps of potential
radiation averaged over time are also computed. The spatial distribution of
time-averaged clear sky solar radiation are calculated over each elevation
band, and *n*_{c} equally distributed radiation classes are identified.
For each radiation class, the mean daily cumulated clear sky radiation value
is computed (termed radiation index, RI, herewith; MJ m^{−2}) and used
in the snowmelt computation. Radiation is pre-processed and provided as model
input. Therefore, the relative module is not included in TOPMELT code, but
made available as a stand-alone tool (see the code availability section at
the end of the paper).

## 2.2 Computation of precipitation amount and phase

Snow accumulation is computed starting from estimates of precipitation and air temperature, based on air temperature and precipitation data from the available weather stations. Similarly to radiation, the model permits the use of several techniques, ranging from Thiessen's polygons to multi-quadratic interpolation (Borga and Vizzaccaro, 1997) for the estimation of basin mean areal precipitation values, which are provided to the model as input data. For the analyses reported in this work, the Thiessen method was used to calculate the mean precipitation over the basin. Air temperature data are used to estimate an unique hourly vertical lapse rate for the whole basin.

To account for gauge catch deficiencies that occur during periods of snow,
precipitation data are corrected with a snow correction factor (SCF). This is
a multiplier of the precipitation data which is applied when station
temperature is lower than a threshold temperature *T*_{c}. Finally, the
basin precipitation value *p*_{basin} is obtained by applying a
non-dimensional precipitation correction factor (PCF) to account for poor
spatial representativeness of rain-gauge stations.

TOPMELT computes the precipitation value at the *i*th elevation band *p*_{i}
(mm h^{−1}) by using a vertical precipitation gradient, accounting for
increased precipitation over elevation. Based on results from Tuo et al. (2016),
this is obtained by means of a precipitation gradient *G* (km^{−1}), as
follows:

where *h*_{i} and *h*_{ref} (m a.s.l.) are the mean altitude of the *i*th elevation band and of the basin respectively. Equation (1) is applied in a way to modify only the distribution of precipitation across the elevation bands without altering the value of *p*_{basin}.

PCF and *G* parameters are generally obtained by comparing model-based snow
cover simulations with satellite-based snow-cover estimates. Since the
average areal precipitation *p*_{basin} is a TOPMELT input, only *G* is
a TOPMELT parameter. Optionally, TOPMELT permits different precipitation
values for each elevation band.

Temperature *T*_{i} (^{∘}C) is provided as input for each elevation band
and time step. In this work, a mean value of air temperature *T*_{i} over the
*i*th elevation band is obtained by using the aforementioned vertical
temperature lapse rate. Estimation of the precipitation phase (solid or
liquid) is therefore performed over each elevation band, according to the
threshold temperature *T*_{c}.

## 2.3 Computation of snow and ice melt

For the generic model cell represented by the *i*th elevation band and the
*j*th radiation class, snowmelt rate *f*_{i,j}(*t*) (mm h^{−1}) at time *t*
is computed taking into account air temperature, clear sky radiation and
albedo. During daytime hours, the snowmelt is given by the following:

where *T*_{i}(*t*) is the elevation band temperature; RI_{i,j}(*t*)
(MJ m^{−2}) is the cell radiation index; *d* (–) is the fraction of
daylight hours in a day (i.e. the number of daylight hours divided by 24);
CMF (mm ^{∘}C^{−1} MJ^{−1} m^{2} h^{−1}) is the combined melt
factor, accounting for both thermal and radiative effects; alb_{i} (*t*) (–)
is the albedo of snow; and *T*_{b}=0 ^{∘}C is a threshold base
temperature. Snow albedo is computed for each elevation band based on
Brock et al. (2000):

where ALBS (–) is the fresh snow albedo, *β*_{2} (–) is a dimensionless
parameter, and ∑_{k}*T*_{i}(*t*_{k}) (^{∘}C) is the sum of the positive
hourly temperatures exceeding the threshold base temperature *T*_{b}
since the last snowfall until the current time *t*.

During nighttime hours, snowmelt is simulated accounting only for air temperature, as follows:

where NMF (mm h^{−1} ^{∘}C^{−1}) is the night melt factor.

For rain-on-snow conditions (Anderson , 1976), melting is computed depending on air temperature and on the energy provided by rain:

where RMF (mm h^{−1} ^{∘}C^{−1}) is the rain melt factor and COST
(^{∘}C^{−1}) is a parameter accounting for the influence of rain on
snowmelt (Carturan et al., 2012). For each model cell, the snow water equivalent
(we_{i,j}; mm) is updated by accounting for snow accumulation,
rain-on-snow, melt and freezing water. Water due to snowmelt or rainfall is
first retained in the snowpack as interstitial water termed liquid water
liqw_{i,j} (mm). When liquid water exceeds a water-holding capacity of the
snowpack (termed LWT), this propagates through the snowpack at a rate DYTIME
(m h^{−1}), to form net water flow at the snowpack base.

When air temperature is less than the threshold base temperature, part of the
liquid water refreezes and liqw_{i,j} is reduced and added to the snowpack
through a freezing rate, termed ice (mm h^{−1}). This is computed as
follows:

where *T*_{b} is the threshold base temperature (Eq. 2) and REFRZ (mm ^{∘}C^{−1} h^{−1}) is the freezing factor.
When we_{i,j} is less than a threshold (termed WETH), ice melt starts. This is computed similarly to snow (Eq. 2), but where the snow albedo is replaced by a constant glacial albedo, ALBG (–), as follows:

At night and during rainfall events, glacier melt is computed by means of Eqs. (4) and (5) respectively.

All model parameters and their values are listed in Table 1,
along with variable names and units. TOPMELT is majorly sensitive to the melt
factor CMF, the most significant calibration parameter of the snowmelt model.
It is constant both in time and space, as the variability of the melt rate is
accounted for by the radiation index and by air temperature. Other important
parameters are the fresh snow albedo, ALBS; the rain melt factor, RMF, which
is a constant calculated from a simplified energy budget (Anderson , 1976);
the precipitation gradient *G*, which drives the re-distribution of
precipitation with the elevation; the base temperature *T*_{b}, which
defines the fusion temperature of snow and ice.

## 2.4 Updating the radiation index distribution: the time-integration routine

As reported in the previous sections, the spatial distribution of clear sky
solar radiation changes with time, based both on astronomic variation of the
radiation flux and its interaction with a complex topography. This implies
that the statistical distribution of the radiation index over each elevation
band will be also modified and should be updated. In general a radiation
class, computed at two different time steps, covers two different areas of
the elevation band. Thus, a pixel belonging to a certain class at a given
time will belong to a different class at another time. Figure 1
shows two different maps of RI, (representing 1 January and 1 April, from an
elevation band ranging from 2000 to 2200 m a.s.l., taken from the basin
selected for the case study of this work. The radiation index is distributed
using 10 equally distributed classes. While the radiation index varies from
1.2 to 20.5 MJ m^{−2} in January, it ranges from 7.9 to
34.7 MJ m^{−2} in April. However, differences are not restricted to the
magnitude of the index. Despite each class having the same area within the
elevation band, their spatial distribution changes from one map to the other.
Examination of the figures shows that a number of pixel belonging to class V
in January are included in class IX in April.

Since the two snowpack state variables, we(*t*) and liqw(*t*), are computed
at the model cell level, pixel transition from a given cell to another must
be accounted for, whenever the radiation index distribution is updated with
time (termed switch date here). To account for pixel transition through
classes, TOPMELT implements an adjusting procedure for model state variables.
The procedure is here described applied to the arbitrary cell state variable
*x*_{i,j}, corresponding to the *i*th elevation band and the *j*th radiation
class of the basin. When updating from one radiation index map to another,
pixels from a certain class can, in principle, move to all other classes, and
pixels from other classes can conversely move to that class. The *x*_{i,j}
variable, which corresponds to certain model cell, should be updated
accordingly. Therefore, a 2-D array accounting for the pixels' transition and
the associated variables among classes is defined, namely the transition
matrix **M**_{i} of the *i*th elevation band.

The transition matrix is *n*_{c}×*n*_{c} sized and is computed
for each elevation band and is unique for each switch date of the radiation
index maps. The element ${M}_{i,j,k}$ of the matrix represents the number of
pixels moving at a switch date from the *j*th source class to the *k*th
destination class, within a given elevation band *i*. ${M}_{i,j,j}$ is the
diagonal element of **M**_{i}, representing the pixels that did not move
from the source radiation class. Provided that the total number of pixels
belonging to a class must remain constant, a property of the transition
matrix is that the sum of the elements along the *i*th row is equal to the
sum of the elements of the *j*th column (i.e. the number of pixels leaving a
class is replaced by an equivalent number of pixels migrating from other
classes):

where *N*_{i,j} is the number of pixels in the *j*th radiation class within
the *i*th elevation band. When TOPMELT switches from one radiation index map
to another, the cell state variable *x*_{i,j} in the new model cell will be
the sum of the pixel contribution from other classes and of the pixels
remaining in the source class, where the number of incoming or remaining
pixels is weighted with respect to the total number of pixels of the source
class. Therefore, *x*_{i,j} is corrected through a matrix **C**_{i} of
correction coefficients, relative to the *i*th elevation band, which can be
derived from **M**_{i} through the following relation:

Following Eqs. (8) and (9) that the sum of each line or
column of the correction factors matrix **C**_{i} must be equal to 1,

The coefficient ${C}_{i,j,k}$ represents the correction factor for the state
variable *x*_{i,j} that must be redistributed among the other classes through
the updating process, within the *i*th elevation band. With the updating, if
*x*_{i,j} is the source class variable band and *x*_{i,j} is transformed
(i.e. destination), the class variable correction is computed through the
following:

or through the equivalent forms:

and

Eq. (13) represents the weighted sum of *x*_{i,k} pixels that
moved from the *n*_{c} source classes to the destination *k*th class.
Since the correction factors matrix **C**_{i} can be computed once for
all simulations, the model
computational efficiency is preserved.

To exemplify the computational flow and its constraints, the example of the
water equivalent (w.e.) state variable is reported here. At a given radiation
index switch, we_{i,j}, will be transferred within the *i*th elevation band
across different classes transforming into ${\widehat{\text{we}}}_{i,j}$, for
$j=\mathrm{1},{n}_{\text{c}}$. The total volume of snow at a given elevation band *i* of
the destination distribution is as follows:

where *A*_{p} is the pixel size. Combining Eqs. (13)
and (14), and provided that the number of pixels is the same for
each class of the *i*th elevation band (${N}_{i,j}={N}_{i}$ for $j=\mathrm{1},{n}_{\text{c}}$),

Eq. (10) plus Eq. (15) yield that the transformed w.e. volume is equal to original volume:

Therefore, Eq. (10) is a constraint that holds conservation of w.e. through the updating process.

## 2.5 Representation of the water equivalent distribution and snow cover

The model allows one to provide the representation of spatially continuous water equivalent maps (as well as any other model cell variable) at a given time. This is carried out by exploiting a routine which links each model cell to the corresponding topographic elements, accounting for variation of the radiation index maps. Then, the water equivalent maps may be converted to snow cover maps by using suitable threshold values. In this work, we used a minimum threshold of 10 mm (Parajka and Blöschl, 2008) for the intercomparison with the MODIS snow cover products.

TOPMELT is integrated within a semi-distributed hydrological model (ICHYMOD; Norbiato et al., 2008), which transforms excess precipitation plus melt contribution into runoff at the outlet of the basin. The total flow routed from TOPMELT to ICHYMOD is the areal weighted sum of each single cell flow, which is made by rainwater and excess meltwater. The model consists of a soil moisture routine and a flow routing routine.

The soil moisture routine uses a probability distribution to describe the spatial variation of water storage capacity across a basin, accordingly to the probability-distributed model (PDM) of Moore (2007). Drainage from the soil enters slow response pathways. The base discharge is routed from groundwater to the catchment outlet through a cubic law storage model. Direct runoff from the proportion of the basin where storage capacity has been exceeded is routed by means of a geomorphology-based distributed unit hydrograph (Da Ros and Borga , 1997), conceptualized by a cascade of two linear reservoirs in series. Runoff from ice melt is transferred to the outlet through two different routes, depending on glacial till imperviousness. Part of the ice meltwater is input to the soil moisture storage, while the remaining fraction flows directly to the outlet as a cascade of two linear reservoirs in series. The base discharge is routed from groundwater to the catchment outlet through a cubic law storage model. Storage-based representations of the fast and slow response pathways yield a spatially lumped representation fast and slow response at the basin outlet which, when summed, gives the total basin flow.

Losses due to evapotranspiration are calculated as a function of potential evapotranspiration and the status of the soil moisture store in the PDM. Potential evapotranspiration is estimated by using the Hargreaves method (Hargreaves and Samani, 1982).

## 4.1 Study site, available data and model set-up

TOPMELT is applied to the Aurino River basin ending at San Giorgio, located in the Adige River system in the eastern Alps, Italy (Fig. 2).

The basin has an area of 614 km^{2}, 2.7 % of which is covered by
glaciers for a total of 16.4 km^{2}. Forest covers 33.5 % of the basin,
pasture and grassland 44.5 %, bare soil and rocks 21.6 % and urban
areas the remaining 0.4 %. Elevation ranges from 817 to 3485 m a.s.l.
Mean annual basin averaged precipitation is around 950 mm, with values
ranging from 850 mm at lower elevations to 1300 mm at the highest
elevations. Precipitation and temperature data at hourly time intervals are
provided by 15 gauging stations (see Fig. 2 for locations), whereas observed
discharge are available at the stream-gauge station in San Giorgio Aurino.
The natural runoff regime is partially altered by the reservoir operations
over the 25 km^{2} Neves Basin. Basin topography is described by means of a
DTM with a 30 m grid resolution.

Satellite observations of snow cover at 250 m resolution have been available for the study basin since January 2011 and are provided by an algorithm based on MODIS observations developed by Notarnicola et al. (2013a, b). With this algorithm, MODIS maps provide the presence of snow, clouds, bare soil or water bodies for each pixel or pixels with no feature detected. 50 MODIS maps are available during the period from 1 January to 30 June 2011 with a percentage of cloud cover of less than 10 %.

The basin was subdivided into 14 elevation bands of 200 m each, ranging from
800 to 3600 m a.s.l. Elevations bands were then subdivided into a number
*n*_{c} of radiation classes. To assess the impact of different spatial
aggregation levels of the radiation index on model results, five types of
class subdivisions of the basin were considered. The elevation bands were
divided into *n*_{c}=1, 5, 10, 15 and 20 classes, yielding five
different spatial aggregation labelled with C1, C5, C10, C15 and C20
respectively. Similarly, to analyse the influence of using aggregation
periods of the radiation index, five different updating times were used for
the computation of the radiation index distribution, with durations of 1, 2,
4, 8 and 12 weeks and labelled W1, W2, W4, W8 and W12 respectively. Variable
temporal and spatial discretization allows for different configurations of
the space–time aggregation of the radiation index. For example, label W4-C10
refers to a model set-up with a temporal aggregation interval of 4 weeks
combined with use of 10 radiation classes per elevation band. It is
interesting to observe that the model set-up W12-C1 resembles a traditional
temperature-index model with a radiation correction for the elevation band
(as in Li and Williams, 2008), whereas the model set-up W1-C20 approximates a fully
spatially distributed implementation of the enhanced temperature-index model.
One should bear in mind that when just one radiation class is used, there is
no need to update the distribution of the snow water equivalent.

## 4.2 The time-integration routine: assessment of pixel transition

An important feature of the model is the use of the time-integration routine to ensure consistency in the snowpack simulation. This routine accounts for pixel transition from one radiation index class to another at the switching time. In this section we analyse the pixel transition by using an index which represents the percentage of migrating pixels over the total number of pixels belonging to a given elevation band, as follows:

where MI_{i} is the migration index of the transition, ${\widehat{N}}_{i}$ is the number of pixels changing class during a switch, and *N*_{i} is the number of pixels of the *i*th elevation band.

The percentage of migrating pixels was computed at four elevation bands: the lowest (from the lowest elevation of the basin) of 817–1000 m; two intermediate bands from 1600 to 1800 m and from 2400 to 2600 m; the highest from 3400 to 3485 m (which is the highest elevation in the basin). The analysis was performed for the five temporal aggregations by using 10 radiation index classes, reporting the mean average migration index over the various switches. Results are reported in Fig. 3a, showing that the percentage of migrating pixels ranges from up to 16 % at W1 temporal aggregation to up to 69 % at W12 temporal aggregation, with a considerable increase of the transition percentage with the increase of temporal aggregation. It is interesting to observe that the transition percentage decreases with increasing elevation of the band, i.e. with decreasing the spatial dispersion of pixels corresponding to a certain class. It is worth noting that the migration index increases with the number of radiation classes (results not shown here for brevity).

Finally, a specific analysis aimed to analyse the magnitude of the transition class change. To highlight this aspect, the percentage of pixels which move by only one class (for example from the second to the third radiation class) was computed and compared to the transition percentage. The percentages were computed and averaged for all the elevation bands. Figure 3b shows this comparison, by considering 10 radiation classes for the five temporal aggregations. For aggregation W12, 42 % of the pixels migrate through one class, 17 % through more than one class and only 41 % do not migrate at all. For aggregation W1, 12 % of the pixels migrate through one class, only 1 % through more than one class and 87 % of the pixels do not move from the source class at the switch time. These results agree with those reported in Fig. 3a and underline the impact of using larger temporal aggregations on the pixel transition between various classes.

## 4.3 Model calibration and validation

The results reported in Sect. 4.2 are obviously independent on the specification of the snowpack model parameters. However, their impact on model results (for instance, on the snow water equivalent spatial distribution) depends on the specification of model parameters.

TOPMELT and ICHYMOD parameters were identified by means of a two-stage procedure, based on comparison of the simulated outflow with the discharge measured at San Giorgio and on comparison of the simulated snow cover with MODIS data, for the period where MODIS data were available. The model set-up W4-C10 was used for the parameter identification. The following statistics were used for comparing simulated and observed discharges:

where *Q*_{obs,t} and *Q*_{sim,t} are the observed and
simulated discharge at time *t*, respectively, ${\overline{Q}}_{\text{obs}}$ is the
average value of the observed discharges and *N* is the number of
observations. Optimal values for bias and Nash–Sutcliffe efficiency (NSE) are 0 and 1, respectively.

To compare the simulated snow cover (SC) area with the MODIS observation, a snow water equivalent threshold of 10 mm was used to declare a snow-covered pixel (Parajka and Blöschl, 2008). Then, the 30 m grids contributing to one MODIS pixel were calculated, and a simulated MODIS-like pixel was considered as snow covered if the percentage of the snow-covered 30 m grid size pixels is equal to or higher than 50 %. MODIS maps with cloud coverage of less than 10 % were used for the analysis. For assessing the correspondence of simulated versus observed values, the accuracy index – ACC – skill measure, based on the contingency table, was used:

where TP are the number of true positives, i.e. where both model and observation agree on the presence of snow on the pixel; TN is the number of true negatives; FN is the number of false negatives, i.e. pixels which are snow covered according to MODIS and where the model simulates no snow; FP is the number of false positives, i.e. pixels which are free of snow according to MODIS and where the model simulates snow. ACC ranges between 0 and 1 with its optimum at 1. Application of a comparison between MODIS data and TOPMELT-simulated snow cover is exemplified in Fig. 4 for a sample date: 6 May 2011. Following Parajka and Blöschl (2008), the accuracy index (Eq. 20) was computed on a pixel base over the 50 cloud-free MODIS maps available. The resulting spatial distribution of the accuracy index is termed the overall accuracy (OA) map.

Model parameter identification was carried out by using data from 1 October 2001 to 30 September 2012 by using the model configuration W4-C10. The period from 1 October 2007 to 30 September 2012 was used for model parameter calibration, with optimization of the statistics bias, NSE and ACC, whereas model validation was carried out over the period 1 October 2001 to 30 September 2007. Model parameter optimization started from the parameterization of the ICHYMOD application by Norbiato et al. (2009). Model error statistics NSE and bias are equal to 0.71 % and 2 %, respectively, for the calibration period, and to 0.71 % and −9 %, respectively, for the validation period.

A graphical comparison of simulated (W4C10) and observed discharges is reported in Fig. 5 for the period May–July 2011, showing the general consistency of the simulation.

The overall accuracy map is reported in Fig. 6b, while Fig. 6b shows the land use over the catchment. The figure shows that simulated snow dynamics agrees (OA > 0.7) with MODIS snow cover detection in 71 % of the area. Lower OA corresponds to forest cover, and the north-facing slopes with forest cover are characterized by very low values of OA. This shows clearly the well-known combined effect of view geometry and forest cover on MODIS snow cover accuracy. Forests make MODIS remote sensing of snow challenging because the presence of trees complicates the monitoring of snow as trees obscure snow on the ground surface (Notarnicola et al., 2013b). View geometry may be a further major error source in MODIS snow-mapping algorithms in forested areas. This is because the gaps in forest canopies, which are essentially the detectable snow fraction in winter, are lower at off-nadir views (Notarnicola et al., 2013b).

## 4.4 Impact of temporal and spatial aggregation on model results

The impact of using different spatial and temporal aggregations on TOPMELT results was carried out by considering both the spatial distribution of the water equivalent and the simulated flow. In this analysis, the finest spatial and temporal discretizations, C20 and W1 respectively, were taken as a reference for the comparisons. For the water equivalent spatial distribution, the assessment was carried out for the period between 1 October 2010 and 30 June 2011, generating a distribution of w.e. at a weekly time step, for a total of 50 simulated snow maps. The intermediate subdivision into radiation classes (C10) was used for the comparison of temporal aggregations W12 to W2 with W1; Additionally, the intermediate temporal aggregation of 4 weeks (W4) was used for the comparison of spatial aggregations C1 to C15 with C20 (i.e. the reference configuration for Fig. 7a is W1-C10, for Fig. 7b is W4C20).

The NSE statistic (Eq. 19) was used to quantify the agreement between the analysis and reference w.e. distribution over space, by excluding from the comparison all the occurrences of zero water equivalence on both maps. One value of the NSE statistic was obtained for each of the 50 maps.

The results are reported in Fig. 7, where the distributions of the NSE statistics are summarized by using box plots. Figure 7a summarizes the results concerning the impact of the temporal aggregation, whereas Fig. 7b shows the analysis concerning the spatial aggregation. Using the longest temporal aggregation (W12) has a noticeable impact on the results, with NSE ranging from 0.88 to 0.99. Even in this case, a considerable improvement is obtained by using a reduced temporal aggregation, such as W4. Interestingly, the results resemble the findings obtained by comparing the model set-up in terms of pixel migration. Using just one radiation class (C1) degrades markedly the accuracy of the w.e. distribution, with NSE ranging from 0.78 to 0.98. Results improve considerably by using five classes (C5).

The impact analysis on flow simulations was carried out by comparing observed and simulated flows over the validation period from 1 October 2001 to 30 September 2007. Results are reported in Table 2 by using the NSE for the various spatial and temporal aggregations. The results are at odds with those reported for the w.e. distribution, showing that the sensitivity of the modelled runoff is negligible. The reported NSE values are very close to each other, ranging from 0.70 to 0.73. The comparison of each runoff simulation to the simulations obtained by using the W4 or C10 reference models, as previously done for the snow water equivalent maps, results in NSE values always in excess of 0.99. The comparison statistics varied only in a limited way when the comparison was focused on the 1 March to 30 June period, i.e. the melting season.

These results are not unexpected. The size of the basin, which largely exceeds the correlation scale of the radiation index, together with the branching nature of the river network, provides indeed a powerful way of averaging out the heterogeneity of snowmelt processes, as shown by Comola et al. (2015), among others.

To better highlight the control exerted by the catchment size on runoff
simulations, we subdivided the study basin into a number of sub-basins
characterized by different drainage areas. We isolated 5 basins with mean
drainage of 20 km^{2}, 10 basins with mean drainage area of 10 km^{2} and
20 basins with mean drainage area of 5 km^{2}. The basins were identified to
ensure the sampling of the whole range of solar radiation characteristics of
the main basin. The model was implemented for these basins by modifying only
those model parameters which depend explicitly on the area and the topography
of the basin, including, of course, the distribution of the radiation index.
The other model parameters were transposed from those identified for the
basin ending at San Giorgio Aurino. In these comparisons, we analysed only
the effect of varying the number of radiation classes. We considered only
three subdivisions for the radiation classes: 1, 10 and 20 classes, and we
used the runoff obtained by using the finest spatial resolution (20 classes)
as a reference against which the other two were contrasted. We used a period
of 4 weeks (W4) as the temporal aggregation. As a matter of fact, the
simulations from the model set-up W4C1 and W4C10 were compared with the
corresponding simulations from the model set-up W4C20. The NSE statistics
were computed only over the 1 March to 30 June period, i.e. the melting
season. The results, reported in Table 3 as an average value of the
NSE statistic over the different basins, show clearly the control exerted by
the catchment size on the effect of using a different number of radiation
classes for runoff simulations. The differences between W4C1 and W4C20 runoff
simulations are considerable (NSE = 0.77) for 5 km^{2} area and rapidly
decrease with an increase of the drainage area. At 10 km^{2} and 20 km^{2}
NSE amounts to 0.91 and 0.99, respectively. On the other hand, no degradation
is reported for runoff simulated by using 10 radiation classes instead of 20,
at least for drainage areas equal to or exceeding 5 km^{2}.

This paper presents TOPMELT, a parsimonious snowpack simulation model which integrates an enhanced temperature-index model into a semi-distributed basin-scale hydrological model. This is obtained by discretizing the full spatial distribution of clear sky potential solar radiation into a number of radiation classes in each elevation band. Snowpack simulation is carried out for each radiation class, rather than for each DTM pixel. This allows one to develop synthesis in modelling approaches to snow simulation and provides the potential for analysing the impact of the spatial and temporal aggregation of radiation fluxes on model results. Furthermore, this approach reduces the computational burden, which is a key potential advantage when parameter sensitivity and uncertainty estimation procedures are carried out.

The impact of temporal and spatial aggregation of the radiation fluxes on
model results is assessed by applying TOPMELT on the 614 km^{2} Aurino River
basin at San Giorgio, in the eastern Italian Alps. The analysis is carried
out by examining five temporal aggregation levels (ranging from 1 to 12
weeks) and five spatial aggregation levels (obtained by subdividing each
elevation band into a number of radiation classes ranging from 1 to 20), with
their impact on the prediction of snow water equivalent distribution and
runoff response.

The assessment of the snow water equivalent simulations clearly shows the degradation of model results when using large temporal and spatial aggregation scales, with model efficiency decreasing up to 20 %. On the other hand, the sensitivity analysis shows that averaging out the radiation index over 4 weeks and using 10 radiation classes as subdivisions has a minimal impact on model results.

Analysis of the runoff response simulations shows that the effects of the
spatial patterns of snow water equivalent are strongly smoothed at the scale
of Aurino Basin at San Giorgio, with minimal deviations over the different
considered model set-up. Examination of TOPMELT-driven runoff results
obtained over internal sub-basins ranging in area from 5 to 20 km^{2}
highlights that the effects of space–time aggregation of the solar radiation
patterns on the runoff response are scale dependent, and that scale
dependency is controlled by the spatial aggregation of the radiation index.
Simulations obtained by using just 1 radiation class degrade in accuracy when
the model is applied over basins of around 5 km^{2}, whereas runoff
simulations obtained with using 10 radiation classes show very limited
degradation over all basin areas considered. These results are important for
driving the implementation of operational applications of TOPMELT. They may
prove to be relevant for the data assimilation of remotely sensed snow cover
information, which may be made more effective with increasing accuracy of
modelled snow cover. They may prove to be relevant also in using the spatial
transferability of enhanced temperature-index model parameters at ungauged
basins.

One important implication of our results concerns the transferability of the
simpler temperature-index model, which is simulated in our cases when TOPMELT
is implemented by using just one radiation class. The results suggest that
this spatial aggregation of the radiation patterns does not impair the
spatial transferability of temperature-index models for runoff simulations of
basins larger than a certain threshold, equal to 5 km^{2} in our case study.

TOPMELT has been developed in Python, version 3.6, and additionally tested with Python 2.7 (Python Software Foundation, https://www.python.org/, last access: 3 December 2019). Python installation requires the following additional modules: datetime, inspect, math, os, sys, pyodbc. The code requires the installation of a SQL database (DB) to store input data and to collects output. TOPMELT was developed and tested with MySQL Community Server (GPL), version 5.7.21. TOPMELT source code and a quick user guide are available at the repository https://doi.org/10.5281/zenodo.1342731 (Zaramella et al., 2018).

MZ was responsible for paper editing, model development and python coding as well as data and simulation analysis. MB carried out paper editing and model development. DZ was responsible for model development. LC carried out the analysis of model application to snow data.

The authors declare that they have no conflict of interest.

The authors would like to acknowledge the time spent on the paper by Massimiliano Zappa and the other two anonymous referees. This paper was improved thanks to their valuable comments and suggestions.

This paper was edited by Bethanna Jackson and reviewed by Massimiliano Zappa and two anonymous referees.

Abudu, S., Sheng, Z., Cui, C., Saydi, M., Sabzi, H.-Z., and King, J. P.: Integration of aspect and slope in snowmelt runoff modeling in a mountain watershed, Water Sci. Eng., 9, 265–273, https://doi.org/10.1016/J.WSE.2016.07.002, 2016. a

Anderson, E. A.: A Point Energy And Mass Balance Model of a Snow Cover. Silver Spring, Md.: U.S. Dept. of Commerce, National Oceanic and Atmospheric Administration, National Weather Service, Office of Hydrology, 1976. a, b

Anslow, F. S., Hostetler, S., Bidlake, W. R., and Clark, P. U.: Distributed energy balance modeling of South Cascade Glacier, Washington and assessment of model uncertainty, J. Geophys. Res., 113, F02019, https://doi.org/10.1029/2007JF000850, 2008. a

Armstrong, R. L., Richard, L., and Brun, E.: Snow and climate: physical processes, surface energy exchange and modeling, Cambridge University Press, ISBN: 9780521 854542, 2008. a

Avanzi, F., De Michele, C., Morin, S., Carmagnola, C. M., Ghezzi, A. and Lejeune, Y.: Model complexity and data requirements in snow hydrology: seeking a balance in practical applications, Hydrol. Process., 30, 2106–2118, https://doi.org/10.1002/hyp.10782, 2016. a

Borga M. and Vizzaccaro A.: On the interpolation of hydrologic variables: formal equivalence of multiquadratic surface fitting and kriging, J. Hydrol., 195, 160–171, https://doi.org/10.1016/S0022-1694(96)03250-7, 1997. a

Brock, B. W., Willis, I. C. and Sharp, M. J.: Measurement and parameterization of albedo variations at Haut Glacier d’Arolla, Switzerland, J. Glaciol., 46, 675–688, https://doi.org/10.3189/172756500781832675, 2000. a

Carenzo, M., Pellicciotti, F., Rimkus, S. and Burlando, P.: Assessing the transferability and robustness of an enhanced temperature-index glacier-melt model, J. Glaciol., 55, 258–274, https://doi.org/10.3189/002214309788608804, 2009. a, b

Carturan, L., Cazorzi, F. and Dalla Fontana, G.: Distributed mass-balance modelling on two neighbouring glaciers in Ortles-Cevedale, Italy, from 2004 to 2009, J. Glaciol., 58, 467–486, https://doi.org/10.3189/2012JoG11J111, 2012. a, b

Cazorzi, F. and Dalla Fontana, G.: Snowmelt modelling by combining air temperature and a distributed radiation index, J. Hydrol., 181, 169–187, https://doi.org/10.1016/0022-1694(95)02913-3, 1996. a, b, c, d

Comola, F., Schaefli, B., Da Ronco, P., Botter, G., Bavay, M., Rinaldo, A., and Lehning, M.: Scale-dependent effects of solar radiation patterns on the snow-dominated hydrologic response, Geophys. Res. Lett., 42, 3895–3902, https://doi.org/10.1002/2015, 2015. a, b

Da Ros, D. and Borga, M.: Use of digital elevation model data for the derivation of the geomorphological instantaneous unit hydrograph, Hydrol. Process., 11, 13–33, https://doi.org/10.1002/(SICI)1099-1085(199701)11:1<13::AID-HYP400>3.0.CO;2-M, 1997. a

Dubayah, R., Dozier, J. and Davis, F. W.: Topographic distribution of clear-sky radiation over the Konza Prairie, Kansas, Water Resour. Res., 26, 679–690, https://doi.org/10.1029/WR026i004p00679, 1990. a

Duffie, J. A. and Beckman, W. A.: Solar engineering of thermal processes, Sol. Energy, 45, 9–17, https://doi.org/10.1016/0038-092X(90)90061-G, 2013.

Essery, R., Morin, S., Lejeune, Y. and B Ménard, C.: A comparison of 1701 snow models using observations from an alpine site, Adv. Water Resour., 55, 131–148, https://doi.org/10.1016/J.ADVWATRES.2012.07.013, 2013. a

Formetta, G., Kampf, S. K., David, O., and Rigon, R.: Snow water equivalent modeling components in NewAge-JGrass, Geosci. Model Dev., 7, 725–736, https://doi.org/10.5194/gmd-7-725-2014, 2014. a

Hargreaves, G. H. and Samani, Z. A.: Estimating potential evapotranspiration, J. Irrig. Drain. E., 108, 223–230, 1982. a

Hock, R.: A distributed temperature-index ice- and snowmelt model including potential direct solar radiation, J. Glaciol., 45, 101–111, https://doi.org/10.1017/S0022143000003087, 1999. a, b, c

Hock, R. and Holmgren, B.: A distributed surface energy-balance model for complex topography and its application to Storglaciären, Sweden, J. Glaciol., 51, 25–36, https://doi.org/10.3189/172756505781829566, 2005. a

Jóhannesson, T., Sigurdsson, O., Laumann, T., and Kennett, M.: Degree-day glacier mass-balance modelling with applications to glaciers in Iceland, Norway and Greenland, J. Glaciol., 41, 345–358, https://doi.org/10.3189/S0022143000016221, 1995. a

Konz, M. and Seibert, J.: On the value of glacier mass balances for hydrological model calibration, J. Hydrol., 385, 238–246, https://doi.org/10.1016/J.JHYDROL.2010.02.025, 2010. a

Klok, E. J., Jasper, K., Roelofsma, K. P., Gurtz, J., and Badoux, A.: Distributed hydrological modelling of a heavily glaciated Alpine river basin, Hydrol. Sci. J., 46, 553–570, https://doi.org/10.1080/02626660109492850, 2001. a

Lee, R.: Forest microclimatology, Columbia University Press, ISBN:9780231041560, 1978. a

Li, X. and Williams, M. W: Snowmelt runoff modelling in an arid mountain watershed, Tarim Basin, China, Hydrol. Proc., 22, 3931–3940, https://doi.org/10.1002/hyp.7098, 2008. a, b

Magnusson, J., Wever, N., Essery, R., Helbig, N., Winstral, A., and Jonas, T.: Evaluating snow models with varying process representations for hydrological applications, Water Resour. Res., 51, 2707–2723, https://doi.org/10.1002/2014WR016498, 2015. a

Moore, R. J.: The PDM rainfall-runoff model, Hydrol. Earth Syst. Sci., 11, 483–499, https://doi.org/10.5194/hess-11-483-2007, 2007. a

Norbiato, D., Borga, M., Degli Esposti, S., Gaume, E., and Anquetin, S.: Flash flood warning based on rainfall thresholds and soil moisture conditions: An assessment for gauged and ungauged basins, J. Hydrol., 362, 274–290, 2008. a

Norbiato, D., Borga, M., Merz, R., Blöschl, G., and Carton, A.: Controls on event runoff coefficients in the eastern Italian Alps, J. Hydrol., 375, 312–325, https://doi.org/10.1016/J.JHYDROL.2009.06.044, 2009. a

Notarnicola, C., Duguay, M., Moelg, N., Schellenberger, T., Tetzlaff, A., Monsorno, R., Costa, A., Steurer, C., and Zebisch, M.: Snow Cover Maps from MODIS Images at 250 m Resolution, Part 1: Algorithm Description, Remote Sens., 5, 110–126, https://doi.org/10.3390/rs5010110, 2013a. a

Notarnicola, C., Duguay, M., Moelg, N., Schellenberger, T., Tetzlaff, A., Monsorno, R., Costa, A., Steurer, C., and Zebisch, M.: Snow Cover Maps from MODIS Images at 250 m Resolution, Part 2: Validation, Remote Sens., 5, 1568–1587, https://doi.org/10.3390/rs5041568, 2013b. a, b, c

Oke, T. R.: Boundary layer climates, Routledge, ISBN:9780415043199, 1992. a

Orgill, J. F. and Hollands, K. G. T.: Correlation equation for hourly diffuse radiation on a horizontal surface, Sol. Energy, 19, 357–359, https://doi.org/10.1016/0038-092X(77)90006-8, 1977.

Östrem, G.: Ice Melting under a Thin Layer of Moraine, and the Existence of Ice Cores in Moraine Ridges, Geogr. Ann., 41, 228–230, https://doi.org/10.2307/4626805, 1959.

Parajka, J. and Blöschl, G.: Spatio-temporal combination of MODIS images – potential for snow cover mapping, Water Resour. Res.,44, W03406, https://doi.org/10.1029/2007WR006204, 2008. a, b, c

Parajka, J., Holko, L., Kostka, Z., and Blöschl, G.: MODIS snow cover mapping accuracy in a small mountain catchment – comparison between open and forest sites, Hydrol. Earth Syst. Sci., 16, 2365–2377, https://doi.org/10.5194/hess-16-2365-2012, 2012.

Pellicciotti, F., Brock, B., Strasser, U., Burlando, P., Funk, M., and Corripio, J.: An enhanced temperature-index glacier melt model including the shortwave radiation balance: development and testing for Haut Glacier d’Arolla, Switzerland, J. Glaciol., 51, 573–587, https://doi.org/10.3189/172756505781829124, 2005. a, b, c, d

Ranzi, R. and Rosso, R.: A physically based approach to modelling distributed snowmelt in a small alpine catchment, IAHS Publ., 205, 141–150, 1991. a

Swift, L. W.: Algorithm for solar radiation on mountain slopes, Water Resour. Res., 12, 108–112, https://doi.org/10.1029/WR012i001p00108, 1976. a

Tuo, Y., Duan, Z., Disse, M., and Chiogna, G.: Evaluation of precipitation input for SWAT modeling in Alpine catchment: A case study in the Adige river basin (Italy), Sci. Total Environ., 573, 66–82, https://doi.org/10.1016/J.SCITOTENV.2016.08.034, 2016. a

Vionnet, V., Brun, E., Morin, S., Boone, A., Faroux, S., Le Moigne, P., Martin, E., and Willemet, J.-M.: The detailed snowpack scheme Crocus and its implementation in SURFEX v7.2, Geosci. Model Dev., 5, 773–791, https://doi.org/10.5194/gmd-5-773-2012, 2012. a

Zappa, M.: Objective quantitative spatial verification of distributed snow cover simulations – an experiment for entire Switzerland, Hydrol. Sci. J., 53, 179–191, https://doi.org/10.1623/hysj.53.1.179, 2008.

Zappa, M. and Pos, F. and Strasser, U., Warmerdam, P., and Gurtz, J.: Seasonal Water Balance of an Alpine Catchment as Evaluated by Different Methods for Spatially Distributed Snowmelt Modelling, Hydrol. Res., 34, 179–202, https://doi.org/10.2166/nh.2003.0003, 2003. a

Zaramella, M., Borga, M., and Zoccatelli, D.: TOPMELT 1.0 (Version 1.0), Zenodo, https://doi.org/10.5281/zenodo.1342731, 2018. a