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.
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.
The equation of the snow heat energy
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
Incoming heat energy
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.
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*.
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%.
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.
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.
The authors are grateful to two anonymous reviewers for their constructive comments and suggestions for improving this paper.