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 *E _{sub}*.

*ɛ*is the deformation of the snow during the time

_{s}*t*, σ

*is the main tension in the snow,*

_{s}*η*is the snow viscosity coefficient and

_{s}*t*is the time.

*is bulk of frozen water, θ*

_{i}*is bulk of water vapor and θ*

_{a}*is bulk of water. Using Equation (2), the average density of snow cover can be defined as: where the average density of snow cover (*

_{w}*ρ*) is determined by the density of frozen water

_{s}*ρ*, the density of water vapor

_{i}*ρ*and the density of liquid water

_{a}*ρ*. All terms of equation are expressed in

_{w}*kg m*

^{−3}.

### The equation of the snow heat energy

*c*is heat capacity at constant pressure,

_{p}*k*is the effective thermal conductivity of snow,

_{e}*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*

^{−1}

*K*

^{−1}), whereas the density and the effective thermal conductivity of snow cover are variable over time.

### Snow density

*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}*SWE*and

*SNO*are highly correlated and that

*ρ*increases according to rising of

_{b}*SNO*. Given that, many authors have shown that the relationship between

*ρ*and

_{b}*SNO*is nonlinear (Pomeroy & Gray 1995; Marchand & Killingtveit 2004; Lundberg

*et al.*2006; Joans

*et al.*2009): where

*a*is the initial density of snow (constant values). Using Equations (5) and (6) the following relation can be defined:

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}.

### Effective thermal conductivity of the snow cover

*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: where

*k*is the effective thermal conductivity and

_{e}*ρ*is the snow density. They reported the fit of Equation (9) with the coefficient of determination (

*R*

^{2}) 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

*Q*is convection,

_{CN}*Q*is condensation,

_{CD}*Q*is radiation,

_{RN}*Q*is heat energy from raindrops, and

_{Rain}*Q*is the total heat energy. The incoming power

*W*is heat energy per unit time

*W*= Σ

*Q*/Δ

_{i}*t*. The incoming heat energy (Σ

*Q*) 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 (

_{i}*λ*). The total melted snow

_{s}*m*(

*kg*/

*ms*) per unit time and unit bulk is: where

*λ*is the energy required to transform 1

_{s}*kg*snow to the liquid state (333,700

*J kg*

^{−1}).

## NUMERICAL METHOD

### Snow accumulation

*SWE*is the snow water equivalent,

*R*daily precipitation sums,

*E*is the daily sublimation rate and

_{sub}*SNO*is the daily snowmelt.

_{melt}*R*is precipitation in the form of rain,

_{R}*T*is the mean daily temperature,

*R*is the daily precipitation sum,

*T*is the bottom limit temperature (−1.1 °C), and

_{S}*T*is the upper limit temperature (3.3 °C). The coefficient of snow occurrence (

_{R}*A*) is defined as:

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.

### Snow melting

*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: 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).

*T*) is equal to the mean daily temperature (

*T*). Mathematical formulation of boundary condition at the snow surface is as follows: where

^{mean}*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

*T*> 0 °C, the Neumann type boundary condition is used (Bartlett & Lehning 2002): 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):

^{mean}*T*= 0 °C which is used to define the temperature above the melting point Δ

_{m}*T*at the contact layer (Δ

*T*=

*T*-

*T*). If Δ

_{m}*T*> 0 °C the snowmelt will occur at a daily rate (Bartlett & Lehning 2002): where

*SNO*is total snowmelt,

_{melt}*c*is the heat capacity of the snow, Δ

_{p}*T*is the temperature above the melting point,

*λ*is the latent heat of melting, and

_{s}*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 (*SNO _{melt}*/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 (

*SNO*). Finally, the snow temperature is obtained utilizing the linear interpolation method: the snow temperature

_{melt}*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 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

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 *Q _{CN}* 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

*Q*was the most intense for the winter season 2004, when air humidity was high and coupled with high temperatures. The share of

_{CD}*Q*in melting process varies along the winter season from 10.9% to 15.9%. The contributing importance of radiation

_{CD}*Q*on the snow melting is in the range 6.8–9.9% and largely depends on the sunshine duration. The heat energy by raindrops

_{RN}*Q*has the less importance in the energy balance. The influence of

_{Rain}*Q*was the highest in the 2010 season when it reached 5.5%.

_{Rain}*R*), sublimation (

*E*) and the

_{SUB}*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).

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.

## 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.