This paper presents one-dimensional numerical model for snowmelt/accumulation simulations, based on the equation of heat energy. It is assumed that the snow column is homogeneous at the current time step; however, its characteristics such as snow density and thermal conductivity are treated as functions of time. The equation of heat energy for snow column is solved using the implicit finite difference method. The incoming energy at the snow surface includes the following parts: conduction, convection, radiation and the raindrop energy. Along with the snow melting process, the model includes a model for snow accumulation. The Euler method for the numerical integration of the balance equation is utilized in the proposed model. The model applicability is demonstrated at the meteorological station Zlatibor, located in the western region of Serbia at 1,028 meters above sea level (m.a.s.l.) Simulation results of snowmelt/accumulation suggest that the proposed model achieved better agreement with observed data in comparison with the temperature index method. The proposed method may be utilized as part of a deterministic hydrological model in order to improve short and long term predictions of possible flood events.

INTRODUCTION

The area of the Earth's surface covered with snow makes up about 14% of the global surface (Lenhing 2005) and the snow cover seems to be a factor in creating the climate (Semadeni-Davies 2003). The role of snow is especially important because it reflects the short-wave solar radiation in a high percentage; for example, in the case of fresh snow, reflection equals nearly 90%. The potential changes of the Earth's snow cover may have significant influence on the global heat energy balance. The reflection of incoming solar radiation represents positive feedback and leads to a reduction of global temperature. The reflected solar energy shall be considerably diminished as a consequence of the snow cover reduction.

From the hydrological point of view, the main characteristic of snow cover is temporal retention capacity of water in the river basin. When the appropriate weather conditions are met, water from the snowmelt increases the runoff. The conditions of water movements depend on the soil moisture and physical characteristics of the soil (Andersen et al. 2001; Anderson & McDonnell 2005; Yusop et al. 2007; Prodanović et al. 2009). Hydrological models often consider the water from the snowmelt as contribution to the surface run-off (Neitsch et al. 2009). In addition, the snowmelt process, accompanied with intense rain events, frequently presents the main cause of fluvial floods (Achleitner et al. 2007).

The snow cover consists of fractions of the air (water vapor) and water, where the water can be represented either by the solid or liquid state. The share of these fractions in the snow cover defines the physical properties of snow. Moreover, it affects physical processes that take place in the snow cover. The heat energy transferred to the snow surface is the most significant driver of snow melting. The distribution of heat energy in the snow is defined by the equation of heat energy, which depends on the density, thermal conductivity and heat capacity of the snow cover.

Several methods for the modeling of processes in the snow cover have been proposed in the literature, the most common ones being: degree-day, temperature index, energy balance and mixed theory. The degree-day method is based on the linear relation between snow melting and the air temperature. This relation can be extended by taking into account the influence of solar energy and can be enhanced using the variable melt rate over time, according to the hydrological season (Stojković & Milivojević 2013; Oreiller et al. 2014). The temperature index method is similar as the degree-day method but uses the sinusoidal changes in melting rate during the hydrological cycle (Melloh 1999); this method presents the constitutive part of several run-off models (SWAT, HBV, SRM, UBC, HYMETI). The most accurate and complex snow melt methods are based on the equation of heat energy (PANTA RHEI, TOPKAPI). These methods are used in many applications which are associated with the snow melting process (Zeinivand & De Smedt 2010; Forster et al. 2014; Shrestha et al. 2015). The well known SAFRAN-CROCUS-IFR model is based on the energy equation and is used for the prediction of avalanches and avalanche warnings in regions of France and Switzerland (Brun et al. 1989). In addition, the methods based on the energy equation are Cold Regions Research and Engineering Laboratory (CRREL) (Jordan 1991), DAISY (Bader & Weilenmann 1992) and SNOWPACK (Bartlett & Lehning 2002). Finally, the mixed theory models are based on the postulate of conservation of ice mass, water and water vapor in natural snow (Morland et al. 1990). Due to the mathematical complexity of this theory for practical application, the snowmelt models are simplified utilizing the following assumptions: there is no snow settling (Gray & Morland 1994), there is no water present in the snow (Gray & Morland 1994) or no metamorphism (Bader & Weilenmann 1992).

The aim of the paper is to propose a relatively simple snowmelt/accumulation model that can be implemented as part of hydrological run-off models, which often consider large areas with limited data. Considering the physical processes in the snow cover, the proposed snowmelt model is based on the equation of heat energy. In order to solve this equation, a novel numerical approach based on variable spatial domain is proposed. This model represents fairly well the complex processes in snow cover such as melting, accumulation and sublimation. The model applicability is demonstrated at the meteorological station (m.s.) Zlatibor in the western Serbia, at 1,028 meters above sea level (m.a.s.l.).

Governing processes in a snow cover

The composition of snowfall varies according to weather conditions, depending, for example, if snowmelt occurs during cold weather with a slight wind or during warmer weather with a snowstorm. The shape of snowflakes depends on temperature and humidity. Transport of snow happens under the windy conditions that contribute to the transport part of the snow cover. Of major importance for the snow cover is solar radiation, which can be divided into short-wave and long-wave radiation. The turbulent exchange between the snow surface and atmosphere produces water movements from the solid state to the vapor. This phenomenon is called the snow sublimation Esub.

Accumulation of snowflakes on the surface leads, subsequently, to their settlement, which takes place in layers so that snowflakes somehow blend into each other. This process of settlement is fast at the beginning of the process, after which it rapidly decreases. The snow is treated as a viscoelastic material with irretrievable deformations (Lenhing 2005). The snow deformation can be expressed by the following equation: 
formula
1
where ɛs is the deformation of the snow during the time t, σs is the main tension in the snow, ηs is the snow viscosity coefficient and t is the time.
The snow cover consists of two phases, water and vapor. The water phase, in accordance with meteorological conditions, consists of ice and liquid water. For fractions of snow cover, the following equation may be written (Lenhing 2005): 
formula
2
where θi is bulk of frozen water, θa is bulk of water vapor and θw is bulk of water. Using Equation (2), the average density of snow cover can be defined as: 
formula
3
where the average density of snow cover (ρs) is determined by the density of frozen water ρi, the density of water vapor ρa and the density of liquid water ρw. All terms of equation are expressed in kg m−3.

The equation of the snow heat energy

The one-dimensional (1-D) differential equation of heat energy is applied to the snow cover profile based on the assumption that the density of snow is homogeneous along the vertical profile. It represents the relationship between snow temperature and incoming heat energy at the snow surface. Mathematical formulation that defines this process is parabolic differential equation (Bartlett & Lehning 2002): 
formula
4
where cp is heat capacity at constant pressure, ke is the effective thermal conductivity of snow, T is the snow temperature, z is the vertical distance measured from the snow top surface, W is the total heat power per unit volume. Heat capacity of snow is a constant value (c = 2,090 J kg−1K−1), whereas the density and the effective thermal conductivity of snow cover are variable over time.

Snow density

The density of fresh snow is between 50 and 200 kg m−3, whereas the density of compact snow is in the range of 300–600 kg m−3 (Mutreja 1986). It is known that density of liquid water is approximately 1,000 kg m−3. The snow water equivalent (SWE) of the snow column is determined by the snow depth (SNO) and its density (ρs): 
formula
5
It is expected that SWE and SNO are highly correlated and that ρb increases according to rising of SNO. Given that, many authors have shown that the relationship between ρb and SNO is nonlinear (Pomeroy & Gray 1995; Marchand & Killingtveit 2004; Lundberg et al. 2006; Joans et al. 2009): 
formula
6
where a is the initial density of snow (constant values). Using Equations (5) and (6) the following relation can be defined: 
formula
7

The initial density of snow as suggested by Pomeroy & Gray (1995) is a = 237 kg m−3. This value is obtained from measurements carried out in North America. Considering the existing concepts and ideas, it is necessary to define a density function from Equation (6) applicable to the climatic characteristics of South-Eastern Europe. The studies that were conducted in Romania gave the result for the mean density and their deviations (Starosolszky 1987). It was shown that the average density of fresh snow varies within the limits 118 ± 28 kg m−3.

For the m.s. Zlatibor, it is assumed that the water equivalent depends on the snow depth according to Equation (7) with the initial snow density as given by the Starosolszky's (1987). The resulting equation is assumed to represent the SWE for the m.s. Zlatibor: 
formula
8

Effective thermal conductivity of the snow cover

The effective thermal conductivity of snow has been frequently studied in Europe and North America. At least 27 studies have been realized with large differences in accuracy in used measurements (Matthew et al. 1997). With data set with 488 measurements, provided by utilization of a thin needle probe inserted in the snow, Matthew et al. (1997) defined the improved regression equations which provide the relationship between effective thermal conductivity and the density of the snow cover: 
formula
9
where ke is the effective thermal conductivity and ρ is the snow density. They reported the fit of Equation (9) with the coefficient of determination (R2) equal to 0.79. It is important to note that the thermal conductivity increases with the snow density. This phenomenon is explained by the fact that water in the solid state has greater thermal conductivity than the air.

Incoming heat energy

The incoming heat energy transferred to the snow surface consists of the following parts (Mutreja 1986): 
formula
10
where QCN is convection, QCD is condensation, QRN is radiation, QRain is heat energy from raindrops, and Q is the total heat energy. The incoming power W is heat energy per unit time W = ΣQit. The incoming heat energy (ΣQi) is transferred to the snow due to an increase of the snow cover temperature, especially at the snow surface layer. As a consequence, the snow melting process occurs. The energy required to transform the snow to a liquid state is the latent heat of melting (λs). The total melted snow m (kg/ms) per unit time and unit bulk is: 
formula
11
where λs is the energy required to transform 1 kg snow to the liquid state (333,700 J kg−1).

NUMERICAL METHOD

Snow accumulation

The snow cover increases when precipitation occurs at negative air temperatures and it decreases with positive temperatures as a result of snow melting and sublimation processes. The balance equation defining the process of change in the SWE over time is: 
formula
12
where SWE is the snow water equivalent, R daily precipitation sums, Esub is the daily sublimation rate and SNOmelt is the daily snowmelt.
Precipitation which occurs near the critical temperature (∼0 °C) is divided into two parts. The first part is retained as snow and the second part is the rain in the liquid state. The distribution of precipitations is obtained with the following linear transformation (Mutreja 1986): 
formula
13
where RR is precipitation in the form of rain, T is the mean daily temperature, R is the daily precipitation sum, TS is the bottom limit temperature (−1.1 °C), and TR is the upper limit temperature (3.3 °C). The coefficient of snow occurrence (A) is defined as: 
formula
14

The range for coefficient A is between 0, for temperature greater than or equal to 3.3 °C, to 1 for temperatures less than or equal to 1.1 °C.

A numerical model of the snow accumulation process is obtained by numerical integration of Equation (12). To solve this equation, the standard Euler method is utilized in the proposed model and the numerical solution for different temperature ranges is given by the following equations: 
formula
15

Snow melting

The process of heat energy transport in the snow column is described by the 1-D differential Equation (4) in the time and spatial domains. The depth of snow cover is changing through the simulation period and, for this reason, a variable spatial domain is considered. The vertical profile of snow cover consists of finite (in this case k = 50) equal lamellas for which the snow temperature is calculated. The physical characteristics of the snow cover along the snow profile are considered as homogeneous, but nonstationary in time. The spatial and time discretization used in this study are adopted as follows: 
formula
16
where SNO is snow depth. Figure 1(a) is an illustration of the discretization of time domain and variable spatial domain in accordance with Equation (16).
Figure 1

Discretization of time domain and variable spatial domain (a); interpolation the snow temperature from time step Δt/2 to Δt (b). T – snow temperature, SNO – snow depth, SNOmelt – snow melting, t – time domain, z – space domain, Δt – time step, Δz – space step; number of lamellas is only for illustrative purposes.

Figure 1

Discretization of time domain and variable spatial domain (a); interpolation the snow temperature from time step Δt/2 to Δt (b). T – snow temperature, SNO – snow depth, SNOmelt – snow melting, t – time domain, z – space domain, Δt – time step, Δz – space step; number of lamellas is only for illustrative purposes.

Solution of the heat energy differential equation requires known initial and boundary conditions. The Dirichlet type boundary condition is utilized at the snow surface, i.e. surface temperature (T) is equal to the mean daily temperature (Tmean). Mathematical formulation of boundary condition at the snow surface is as follows: 
formula
17
where t and z present time and space domain and z = min corresponds to the snow surface at the time step. This boundary condition is applied if the mean daily temperature is less than or equal to 0 °C. If Tmean > 0 °C, the Neumann type boundary condition is used (Bartlett & Lehning 2002): 
formula
18
The Dirichlet type boundary condition is applied at the surface layer between the snow cover and the Earth. The temperatures of the Earth's surface do not differ much from the 0 °C (Bartlett & Lehning 2002): 
formula
19
It should be noted that snow temperature cannot be greater than 0 °C. The excess of heat energy initiates the snowmelt and after that snow temperature will be adjusted to 0 °C at the contact layer. The melting temperature is Tm = 0 °C which is used to define the temperature above the melting point ΔT at the contact layer (ΔT = T-Tm). If ΔT > 0 °C the snowmelt will occur at a daily rate (Bartlett & Lehning 2002): 
formula
20
where SNOmelt is total snowmelt, cp is the heat capacity of the snow, ΔT is the temperature above the melting point, λs is the latent heat of melting, and m is the mass of snow in the vertical profile at the time step t.

In order to present highly dynamic natural process such as transport of the heat energy along the snow column, the snow temperature is calculated at the half time interval (Δt/2). The reason for using this approach is that the daily amount of the snowmelt rate can reach significant values and, therefore, it is important to be estimated with higher accuracy. The hour temperature oscillations are not taken into account for calculation of the snow temperature. Interpolation of the snow temperature at the half time interval (Δt/2) towards the end of the time interval (Δt) is conducted in several steps, as illustrated in the Figure 1(b). In the first step, the snow depth SNO at time step t=m is reduced by half melted snow during the time t = m + 1 (SNOmelt/2). Then, the snow depth SNO and temperature T are calculated at the time step t = m + 1/2. In the next step, the snow depth SNO (t=m + 1) is determined using the whole rate of the snowmelt (SNOmelt). Finally, the snow temperature is obtained utilizing the linear interpolation method: the snow temperature T from the time step t = m + 1/2 is spatially interpolated to the time step t = m + 1. The interpolated temperature is denoted as T*.

The 1-D differential equation of heat energy in the snow cover is solved using the implicit approach. The finite difference numerical formulation of Equation (4) is given the following form: 
formula
21

The implicit parameter Θ in Equation (21) is taken as 0.75. The system of algebraic equations provides a solution to determine the snow temperature profile at specific time step. The system of linear equations is three-diagonal matrix which is solved by the Thomas algorithm (Hoffman 2001).

CASE STUDY

The proposed model is applied for simulation of the snowmelt/accumulation at m.s. Zlatibor, situated in the western mountainous region of Serbia (Figure 2). Spatial coordinates of the location are 43°44′ and 19°43′, with elevation 1,028 m.a.s.l. The meteorological station is part of the network system of the Hydro-Meteorological Service of Serbia (HMSS). HMSS performs a continuous, official monitoring of weather conditions through measurements and observations of numerous parameters.
Figure 2

Location of m.s. Zlatibor (43°44′, 19°43′, 1,028 m.a.s.l.) with the hydrographic network of Serbia.

Figure 2

Location of m.s. Zlatibor (43°44′, 19°43′, 1,028 m.a.s.l.) with the hydrographic network of Serbia.

For simulation of the snow melting process, daily record series being used include air temperature, precipitation and snow depth. Furthermore, the additional daily meteorological parameters such as extreme temperatures, vapor pressure, air humidity, solar radiation, isolation and wind intensity are utilized to determine incoming heat energy by using Equation (10). The heat energy includes convection, condensation, radiation, and heat energy from raindrops. By utilizing the estimated heat energy on the snow cover, the Neumann type boundary condition is applied to calculate the snow temperature and snow melt rate by Equations (18) and (20), respectively. Consistently, the snow melt and sublimation rates are subtracted from the mass of the snow profile, according to Equation (15). At the final step, the 1-D differential equation of heat energy (Equation (21)) is solved in order to determine the distribution of temperature in the snow profile.

RESULTS AND DISCUSSION

The total heat energy which is transferred to the snow surface is calculated by Equation (10). The rapid melting of the snow cover is mainly caused by heat energy from convection QCN associated with a significant increase in the temperature. The contributing importance of convection in the melting process is in the range of 72.1–78.1%. Decreasing snow cover by condensation QCD was the most intense for the winter season 2004, when air humidity was high and coupled with high temperatures. The share of QCD in melting process varies along the winter season from 10.9% to 15.9%. The contributing importance of radiation QRN on the snow melting is in the range 6.8–9.9% and largely depends on the sunshine duration. The heat energy by raindrops QRain has the less importance in the energy balance. The influence of QRain was the highest in the 2010 season when it reached 5.5%.

Figure 3(a) shows the cumulative values of daily precipitation (R), sublimation (ESUB) and the SWE for m.s. Zlatibor for four winter seasons. The SWE is calculated using the proposed numerical model of snow accumulation and melting according Equations (15) and (21).
Figure 3

Cumulative precipitation R, sublimation Esub and SWE (a); the simulated SWE based on the 1-D equation of the heat snow energy (EB) and the temperature index method (TI) for four winter seasons at m.s. Zlatibor (b).

Figure 3

Cumulative precipitation R, sublimation Esub and SWE (a); the simulated SWE based on the 1-D equation of the heat snow energy (EB) and the temperature index method (TI) for four winter seasons at m.s. Zlatibor (b).

The processes in the snow cover are modeled by using the energy balance equation. Furthermore, the proposed method is compared with the temperature index method which is a part of many run-off hydrological models (Stojković & Milivojević 2003; Prodanović et al. 2009). The results are shown in Figure 3(b).

As a parameter of the model performance the correlation coefficient (r) is utilized. The values of r are in the range of 0.906–0.947 and 0.750–0.886 for all considered seasons for the energy balance method and the temperature index method, respectively (Figure 3(b)). It can be concluded that the energy balance method has achieved significantly better agreement with observed values of the snow cover depth. Hence, it is reasonable to assume that the proposed methodology is more appropriate for utilization as a component of hydrological run-off models.

The temperature index method shows higher values of accumulated snow cover than records during the four seasons (Figure 3(b)). Also, the melting process is shifted towards the last time period in the spring and it occurred more rapidly than the method based on the equation of heat energy. This phenomenon is a consequence of meteorological factors such as wind, solar radiation, humidity, and raindrops not being considered.

When the temperature is greater than 0 °C, the excess of heat energy is used to melt the snow cover. Given that, the melting rate reduces the snow cover. The proposed model calculates the vertical distribution of the snow temperature by using the Crank–Nicolson method. Figure 4 illustrates the snow temperature at the vertical profile of the snow cover for the winter season 2003 at the m.s. Zlatibor for several time sections.
Figure 4

Snow temperature distribution T (°C) at the vertical profile of the snow cover (SWE) obtained by using the Crank–Nicolson method at m.s. Zlatibor for the selected time steps t = 9, 26, 29, 34, 59, 69 for the winter season 2003.

Figure 4

Snow temperature distribution T (°C) at the vertical profile of the snow cover (SWE) obtained by using the Crank–Nicolson method at m.s. Zlatibor for the selected time steps t = 9, 26, 29, 34, 59, 69 for the winter season 2003.

CONCLUSIONS

In this study, a novel numerical method for snow accumulation/melting has been developed based on 1-D heat energy. The space domain changes along the time steps and the snow column is discretized by the constant number of numerical nodes. Snow characteristics, such as density and thermal conductivity, are assumed constant over the snow depth, but vary with time. The numerical solution of the heat equation is obtained by using the finite difference implicit numerical scheme. The application of the proposed model is demonstrated by simulation of the snow cover observed at the m.s. Zlatibor in Serbia during several cold seasons. The results have shown very good agreement with the observed snow depths. In addition, the calculated results have been compared with the widely utilized temperature index method. The comparison indicates that the deviation from records is considerably greater when using the temperature index method. Based on these results, and considering the model simplicity, it is reasonable to expect that the proposed numerical method is appropriate for utilization in hydrological run-off modeling, in the case of limited data availability.

The surface run-off generation by the snowmelt is mainly controlled by melting and refreezing processes of the snow cover. Refreezing processes are most significant at the surface layer, especially in the cold regions (Pelt et al. 2012). Hence, the proposed numerical method can be modified to model this process by using an hourly time scale instead of a daily one. In this case, inter-daily temperature variability can have significant impact on the model accuracy.

ACKNOWLEDGEMENTS

The authors are grateful to two anonymous reviewers for their constructive comments and suggestions for improving this paper.

REFERENCES

REFERENCES
Andersen
H. E.
Pedersen
M. L.
Jorgensen
O.
Kronvang
B.
2001
Analysis of the hydrology and flow of nitrogen in 17 Danish catchments
.
Water Sci. Technol.
44
(
7
),
63
68
.
Anderson
M.
McDonnell
J.
2005
Encyclopedia of Hydrological Sciences: Rainfall-Runoff Modeling
.
The Aritrium, Sautern Gate
,
Chichester
,
UK
.
Brun
E.
Martin
E.
Simon
V.
Gendre
C.
Coleou
C.
1989
An energy and mass model of snow cover suitable for operational avalanche forecasting
.
J. Glaciol
.
35
(
121
),
333
342
.
Gray
J. M.
Morland
L.
1994
A dry snow pack model
.
Cold Reg. Sci. Technol.
22
(
2
),
135
148
.
Hoffman
J.
2001
Numerical Methods for Engineers and Scientists
, 2nd edn.
Marcel Dekker
,
Inc
.,
New York, NY
,
USA
.
Jordan
R.
1991
A One Dimensional Temperature Model for a Snowcover
,
CRREL Special Report. 91–16
.
US Army Corps of Engineers
,
USA
.
Lenhing
M.
2005
Energy Balance and Thermophysical Processes in Snowpacks.
Encyclopedia of Hydrological Sciences: Rainfall-Runoff modeling
.
Wiley, Chichester
,
UK
.
Lundberg
A.
Richardson Naslund
C.
Andersson
C.
2006
Snow density variations: consequences for ground-penetrating radar
.
Hydrological Processes
20
(
7
),
1483
1495
.
Marchand
W. D.
Killingtveit
A.
2004
Statistical properties of spatial snowcover in mountainous catchments in Norway
.
Nordic Hydrology
35
(
2
),
101
117
.
Matthew
S.
Holmgren
J.
Konig
M.
1997
The thermal conductivity of seasonal snow
.
J. Glaciol.
43
(
143
),
26
41
.
Melloh
A. R.
1999
A Synopsis and Comparison of Selected Snowmelt Algorithms
.
CRREL Report
.
US Army Corps of Engineers
,
USA
.
Morland
L. W.
Kelly
R. J.
Morris
E. M.
1990
A mixture theory for a phase-changing snowpack
.
Cold Reg. Sci. Technol
.
17
(
3
),
271
285
.
Mutreja
K. N.
1986
Applied Hydrology
.
Tata McGraw-Hill
,
New Delhi
,
India
.
Neitsch
S.
Arnold
J.
Kiniry
J.
Williams
J.
2009
SWAT: Solar and water assessment tool 2009–Theoretical Documentation
.
Texas Water Resource Institute Technical Report No. 406
.
Texas A & M University System, College Station, Texas
,
USA
.
Oreiller
M.
Nadeau
D. F.
Minville
M.
Rousseau
A. N.
2014
Modelling snow water equivalent and spring runoff in a boreal watershed, James Bay, Canada
.
Hydrological Processes
28
(
25
),
5991
6005
.
Pelt
W. J. J.
Oerlemans
J.
Reijmer
C. H.
Pohjola
V. A.
Pettersson
R.
Angelen
J. H.
2012
Simulating melt, runoff and refreezing on Nordenskioldbreen, Svalbard, using a coupled snow and energy balance model
.
Cryosphere
6
,
641
659
.
Pomeroy
J. W.
Gray
D. M.
1995
Snow Accumulation, Relocation and Management
.
National Hydrology Research Institute Science Report No. 7
.
Saskatoon
,
Canada
.
Prodanović
D.
Stanić
M.
Milivojević
N.
Simić
Z.
Stojanović
B.
2009
Modified rainfall-runoff model for bifurcations caused by channels embedded in catchments
.
J. Serbian Soc. Computat. Mech.
3
(
1
),
111
126
.
Semadeni-Davies
A.
2003
Response surfaces for climate change impact assessments in urban areas
.
Water Sci. Technol.
48
(
9
),
165
175
.
Starosolszky
Ö.
1987
Applied Surface Hydrology
.
Water Resources Publication
,
Littleton, CO, USA
.
Stojković
M.
Milivojević
N.
2013
Hydrological modeling with special reference to snow cover processes
.
Facta universitatis–series Architecture and Civil Eng.
11
(
2
),
147
168
.