1 Introduction
The amount of snow, or snow depth, is important to many applications like road safety (Juga et al.; 2013), risk assessments (Bocchiola et al.; 2006; Blanchet and Davison; 2011), winter sport activities (Snow Forecast; n.d.; Fnugg; n.d.), hydrology (Jonas et al.; 2009) and climate change research (Brown et al.; 2003; McCabe and Wolock; 2010; Falarz; 2002; Scherrer et al.; 2013; Zhang et al.; 2004).
The modeling of snow depth is typically divided into three parts: snow accumulation (snowfall), snow aging and melting. The physics behind the different parts is quite complicated and depends on many factors. For snowfall, a common rule is the 10:1 rule stating that the density of the arriving snow is one tenth of the density of water. Following this rule, 10 mm of precipitation results in 10 cm of snow. In reality, the relation is more complicated. The density of snowfall is related to the icecrystal structure by virtue of the relative proportion of the occupied volume of crystal composed of air. Snow density is regulated by incloud processes that affect the shape and size of growing ice crystals, subcloud processes that modify the ice crystal as it falls, and groundlevel compaction due to prevailing weather conditions and snowpack metamorphism. Understanding how these processes affect snow density is difficult because direct observations of cloud microphysical processes, thermodynamic profiles, and surface measurements are often unavailable. Roebber et al. (2003)
builds a neural network to classify snow density of snowfall to three classes heavy (1:1
ratio 9:1), average (9:1 ratio 15:1), and light (ratio 15:1) where ratio refers to the density of water compared to the density of the arriving snow. The authors use the predictors solar radiation, temperature, humidity, precipitation and wind. The method only classify to correct density class in 60% of the cases which emphasize the difficulty in forecasting snow density and snow depth. Snow aging and melting of snow and resulting changes in snow depth depends on many factors like temperature, precipitation, solar radiation and the age and density of the snow pack. Compared to other factors, precipitation and temperature are the main drivers of changes in snow depth and most snow models are only based on these factors (Brown et al.; 2003; Kohler et al.; 2006).Weather forecast services routinely forecast quantities like temperature, precipitation, wind and air pressure, but very rarely snow depth or other snow related quantities. An exception is snowforecast.com
(Snow Forecast; n.d.) which forecast snow depths for skiing resorts around the world, but the methods used are not publicly available.
The small number of weather services forecasting snow depth or other snow related quantities are in a big contrast to the amount logged snow depth data available around the world. E.g. in the US and Canada daily snow depth are logged for at least 8000 locations (Brown et al.; 2003). Motivated by the lack of models to forecast future snow depths and the large amount of historic snow depth data available, in this paper we present models to fill this gap. We build both shortterm forecasts, where reliable forecasts of temperature and precipitation are available, and forecasts going further into the future. The work in Brown et al. (2003) is related to the work in this paper, except that the authors use a numerical model to relate snow depth to precipitation and temperature, while we present a statistical model such that probabilistic forecasts can be performed.
2 Short term forecasting of snow depth
Temperature and precipitation are the main drivers of changes in accumulated snow depth. Better forecasts of future snow depth can therefore be achieved by including weather forecasts of temperature and precipitation in the forecast model. Such a model will be presented in this section.
Let denote the current snow depth at a specific time of day at day where is the number of days with observations. Further let and denote the total precipitation and average air temperature for the last 24 hours.
Snow depth will for a large portion of days be zero (no snow), and else always larger than zero. Precipitation has the same properties and is typically modeled by a zeroinflated gamma model (Stern and Coe; 1984; Sloughter et al.; 2007; Möller et al.; 2013). The model results in a good fit also to the snow depth data and can be formulated as follows
(1) 
where is the gamma density and return if is true and 0 else.
We model as a generalized linear model, and the expectation of the gamma density is linked to covariates as follows. We follow the typical assumption of changes in snow depth by dividing in snowfall (increase) and snow aging/snow melt (decrease) (Kohler et al.; 2006; Kuusisto; 1984; Brown et al.; 2003). First we find a model for snowfall. If the temperature is sufficiently cold, precipitation will arrive as snow, and if the temperature is sufficiently warm, the precipitation will arrive as rain. For temperatures around zero C, precipitation will arrive as a mixture of snow and rain. Several functions are suggested to model this, see e.g. Wen et al. (2013); Kienzle (2008)
. We use the inverse logit function which fits well to our data and is also in accordance with observations from earlier research, see e.g. Figure 6 in
Kienzle (2008). More specifically we assume that the expected depth of snow from mm of precipitation at temperature C is given by(2) 
where
(3) 
The parameter refers to the snow water density ratio and is normally set to 10 (Roebber et al.; 2003)
. We will estimate the parameter from snow depth observation. We expect that
will be negative such that higher temperatures result less snow (lower snow depth).The snowpack tend to sink with time, but less if it is very cold. Further, for temperatures around and above zero C the snow will additionally melt (transformed to water). The aging/melting is going faster if it is raining on the snowpack and is typically modeled with the index method (see e.g. Scherrer et al. (2013)) which simply is an interaction between the temperature and the amount of rain, . To ensure nonnegative snow depths we insert the index method in an inverse logit function. This means that we assume that the expected portion of the snow (depth) that disappears is given by
(4) 
Adding the two parts together, given that we assume that the expected snow depth at time , is given by
(5) 
Since the expectation in (5
) always is larger than zero, we simply use the identity link function when linking the expectation to the gamma distribution. The intercept
typically becomes very small (see Table 2).In gamma regression the far most common is to assume that the shape parameter, , is constant such that the coefficient of variation becomes constant
(6) 
where
denotes the standard deviation of
. Under this parameterization the standard deviation of increases with the expected value of . This is not a natural assumption for the modeling snow depth and turns out to give a very poor fit to the data. If the temperature is below zero C and it is no precipitation, the changes in snow depth is small, and can be predicted with high precision from even whenis high. A more natural assumption is that the variance depends on the expected
change in snow depth. The larger expected changes, the larger uncertainty. To avoid the possibility that the variance is equal to zero, we combined this with the assumption of constant variance resulting in the following model for the variance(7) 
Given the expectation and variance, the shape, , and scale, , parameter of the gamma distribution can be computed in the usual way
(8) 
An alternative to the model formulation above is to fit the data with a double generalized linear model, see e.g. the dglm
package (Dunn and Smyth; 2015) in R (R Core Team; 2015). We have not tried this.
Next we turn to in (1
) which is modeled by logistic regression. One may estimate this probability using the covariates
, and . Instead we use as the only covariate(9) 
which turns out to perform well. The intuitive is that higher expected snow depth, given by (5), results in lower probability of no snow.
Finally, given the covariates and for , we assume that and are independent. This makes it straight forward to put up the likelihood function for the snow depth observations. Due to the nonlinearities to the covariates and the coupling between the logistic and gamma part of the model, to the best of our knowledge there exists no statistical packages in R (R Core Team; 2015) or elsewhere to estimate the parameters of the model. Instead we estimate the parameters by implementing a steepest descent optimization algorithm and find the parameters the optimize the likelihood function.
Given the estimated parameters of the model, forecasting of future snow depths can be computed using Monte Carlo simulations. We assume that a reliable weather forecast of and
is available for the next few days. The model above can then be used to “track” the probability distribution of snow depths into the future as follows. Given the current snow depth
and weather forecasts of temperature and precipitation the next day ( and ), generate a large set of realizations from the distribution in (1). Next, for each sample from and given weather forecasts two days into the future ( and ), generate a sample from . In such a way we track the snow depth into the future given the weather forecasts of temperature, , and precipitation, (.We tried some extensions to the model as described below, but neither of them improved the model with respect to the Akaike information criterion (AIC).

One may expect an unsymmetry for the snowfall inverse logit function and thus we considered the extension .

For the melting part of the model one may argue that if is large, not the whole snow pack will be exposed for the air temperature and snow melting may go slower. We therefore tested the extension of the melting part of the model .

We also conditioned on previous values of snow depth adding the term to (5).
3 Long term forecasting of snow depth
Weather forecasts for temperature and precipitation are typically reliable for three to six days into the future. In this section, we build models to forecast snow depth further into the future than this timespan. We consider two different strategies

Model 1: When reliable weather forecasts are not available, we use historical observations of precipitation and temperature to build statistical time series models for these variables. The forecasts from these models are further used as input to the model in the previous Section.

Model 2: We build a time series model for the snow depth data directly.
The strength of model 1 is that the model gives us simultaneous forecasts of temperature, precipitation and snow depth trends.
3.1 Model 1
Let denote the day during a season for observation time . E.g. if refers to December 31 for some year, for ordinary years and 366 for leap years. For longterm forecasts of temperature and precipitation, the seasonal trends will be important. We apply Fourier series, which are able to model complex seasonal patterns with only a few parameters
(10) 
Except for the seasonal trend, temperature data fits well to an autoregressive process. Thus, we model the temperature time series using an Autoregressive process with a seasonal trend given by (10)
(11) 
where where
denote a normal distribution with expectation
and standard deviation . Further we assume that are independent. Standard packages in R Core Team (2015) can be used to fit this model, but instead we found the parameters that maximized the likelihood function using a steepest descent optimization algorithm.For precipitation we follow the model in Stern and Coe (1984) with some exceptions that will be explained below. Because of the many days with no precipitation the zeroinflated gamma model in (1) is suitable also to model precipitation
(12) 
Similar to the model for snow depth and the model in Stern and Coe (1984) we model with a gamma regression model. We define the expectation as
(13) 
were is defined such that if and if . Stern and Coe (1984) also considers interactions between for different values of . We achieve almost as good fit with respect to AIC by instead increasing the value of . The difference in AIC using interactions or not were between 2 and 8 for the data series considered in this paper. For simplicity, we therefore omitted interactions, which made the model easier to interpret. In contrast to Stern and Coe (1984) we also include the current temperature as a predictor, which results in a substantial improvement of the model. The standard approach of holding the coefficient of variation constant (equation (6)) resulted in a good fit to the precipitation data.
In Sloughter et al. (2007); Möller et al. (2013) it is suggested to model in stead of with a gamma distribution. Based on goodness of fit analyzes we were not able to show that one of these alternatives resulted in a better fit then the other and decided to model as gamma distributed as shown above.
Similar to the snow depth model, we model with logistic regression. For the snow depth data using as the only covariate performed well, but using only as a covariate turns out to perform poorly for the precipitation data. A better fit is achieved using the same covariates as above
(14) 
Also for this model, we used a steepest descent optimization algorithm to find the parameters that maximized the likelihood function.
Given forecasts of temperature and precipitation using the models above, the model in Section 2 can further be used to forecast snow depth.
3.2 Model 2
While model 1 in the previous section perform long term forecasts of snow depth by first doing long term forecasts of temperature and precipitation, in this section we instead model the snow depth time series directly. The model will be exactly as the model in Section 2 except that the covariates must be changed since precipitation and temperature is unknown. Therefore we change (5) with
(15) 
were is defined such that if and if . It turned out that using both the covariates and resulted in a better fit than using only or only .
3.3 Monte Carlo procedure
For the first days in to the future when reliable weather forecast of temperature and precipitation is available, the Monte Carlo method described in Section 2 will be used. Let denote the number of days with reliable weather forecasts. Now suppose that we want to forecast . After running the Monte Carlo method in Section 2 we have a set of realizations of the time series . Longterm forecasts can then be computed as follows

Using model 2, forecasting of can be achieved by generating a realization from model 2 conditioned on each of the realizations of the time series .

Using model 1, first a large set of realizations of and is generated conditioned on the weather forecasts and . For each realization of and , a realization of is generated using the model in Section 2.
The procedures above can be repeated for as long into the future that forecasts of snow depth is needed.
3.4 Goodness of fit
Assume that
is a stochastic variable with a cumulative distribution function
. It’s a wellknown fact that wheredenote a uniform distribution on the
interval. The procedure can be used for the different models presented above. For the model in Section 2 the cumulative distribution for can be computed as(16) 
where is the cumulative gamma distribution. Let denote the real observations of snow depth. If the model fits the data well, we expect the distribution of to be close to uniformly distributed. In the computation of , the real observations of temperature, precipitation and snow depth from the previous day are used as input according to (5). The same procedure will be used also for the other models above.
4 Real data example
In this section, we forecast future snow depths using the models introduced in the previous sections. We downloaded average daily temperature, precipitation and snow depth from three locations in Norway from the web portal eklima.met.no
. Properties of the three locations are shown in Table 2.
Place  Altitude (m)  Timespan with observations 

Oslo (Blindern)  97  June 11 1955 June 11 2015 
Geilo  810  September 1 1966 November 30 2006 
Tromsø  100  January 1 1955 September 1 2015 
The locations represent very different climates in Norway. Oslo and Tromsø are close to the coast in the southeast and far north of Norway respectively, while Geilo is far inland in the mountainous parts of Norway.
We now fit the model in Section 2 for each of the three locations using the time series of temperature, precipitation and snow depth. Table 2 shows the estimated parameters for the three locations.
Place  

Oslo  0.96  0.88  1.76  1.99  0.30  0.03  4.13  1.97  0.63  1.79  
Geilo  0.72  1.74  1.19  2.86  0.25  0.05  3.61  0.75  1.04  2.83  
Tromsø  0.89  2.02  0.83  2.64  0.16  0.04  3.38  0.64  0.98  8.56 
For the fitted models Figures 1, 2 and 3 show curves for for different values of , and . are computed using the law of total expectation
(17) 
Figure 1 snows expected snow depth from 10 mm precipitation for different temperatures. We set so that we only look at the snowfall part of the model, recall (5).
As expected the snow depth decreases rapidly around C and the curves are in accordance with earlier research on such curves (Wen et al.; 2013; Kienzle; 2008). Due to differences in climate the curves varies with location in Norway and such differences is also observed in earlier research (Wen et al.; 2013; Kienzle; 2008).
Figure 2 shows for different temperatures when assuming that mm and cm. Since is set to zero, this figure shows the melting/aging part of the model.
Figure 3 shows for different amounts of precipitation when C and cm. This figure shows the effect of precipitation on the reduction of snow depth. We see that rain has a strong impact on the reduction of snow depth. We also observe that the curves seem to be quite linear which is in accordance with the index model (Scherrer et al.; 2013).
The Figures 1 3 show that the model captures the main effects of snow accumulation, aging and melting.
To perform long therm forecasts of snow depth, we fitted the time series of temperature, precipitation and snow depth using the models in Section 3. In each of these models the number of covariates where chosen based on AIC in a forward stepwise regression procedure. E.g. for temperature, first is increased to one, then to one, then to two and so on. The results are shown i Table 3.
Temperature  Precipitation  Snow depth  
Place  
Oslo  2  3  3  5  4  3  5  4  3  1  5 
Geilo  2  3  2  3  3  3  5  3  3  1  3 
Tromsø  2  4  3  6  3  2  6  3  3  1  6 
4.1 Goodness of fit
Figure 4 shows goodness of fit histograms of the models in Sections 2 and 3 using the approach described in Section 3.4. The upper row in Figure 4 shows from left to right goodness of fit histograms for temperature and precipitation, respectively, while the bottom row shows from left to right goodness of fit histograms for the model in Section 2 and model 2, respectively. The goodness of fit histograms for temperature and precipitation are based on observation for the whole year while for the histograms for the snow depth models are based on the winter months December, January and February. All the histograms are for Oslo, but the histograms for Geilo and Tromsø were similar.
We see that the histograms for temperature and precipitation looks fairly uniformly distributed. The two models for snow depth (bottom row) show an overrepresentation of values around 0.5. When the temperature is below C and it is no precipitation, the changes in snow depth is minimal and the model overestimates the variance in this cases resulting in to many values around 0.5 in the goodness of fit histograms. During winter, these weather conditions is of course very common. We can of course reduce the value of in (7) to reduce the variance when the expected change in snow depth is small, but that will result in other negative consequences for the model.
4.2 Forecast performance
We now inspect the forecast performance of the snow depth models. Forecasts are performed using the Monte Carlo procedures described in Sections 2 and 3. In each time step we forecast using the mean value of the Monte Carlo samples. Classification error is measured by the average difference in absolute value between the observed and forecasted snow depth. Classification is performed for the winter months December, January and February in a cross validation procedure. All data except for one year from July 1 to June 30 the next year is used to fit the models, and further used to forecast snow depths for December, January and February for the year of data not included in the model fitting. The procedure is repeated for each year. We consider three different cases where we assume that reliable weather forecasts of temperature and precipitation are available for zero, five and ten days into the future. For the days were we assume that reliable weather forecasts are available, the real observed values of temperature and precipitation are used, i.e. we assume “perfect” weather forecasts. For comparison we also consider a version of model 1 and 2, where all covariates except the periodic covariates are set to zero, i.e. all the variables , , , , , , are set to zero which means that only , , and can larger than zero. In addition, we assume that the number of days with reliable weather forecasts of precipitation and temperature are zero. In other words, these versions of model 1 and 2 do not take advantage of previous observations for the given season or weather forecasts and forecast only based on seasonal properties. We expect that the further we forecast into the future, the less will the usefulness of previous observations for the given season or weather forecasts be.
The gray and black curves show results for model 1 and model 2, respectively. The dashed, dotted and dashdotted curves show forecasts when we assume that reliable weather forecasts of temperature and precipitation are available for zero, five and ten days, respectively. For the days when reliable weather forecasts are available, the model in Section 2 is used. The solid curves show classification performance where only the periodic covariates are included as explained above. The left vertical axis show average classification error in absolute value, while the right vertical axis show the classification error normalized with the average snow depth for the months of December, January and February.
We see that model 2 perform better then model 1. Further we see that given reliable weather forecasts the classification error is about half compared to not having reliable weather forecasts. We also observe that useful forecasts is possible long in to the future. Forecasting three weeks into the future and given five days of reliable weather forecasts the classification error is a little over half of the classification error using only the periodic covariates. In comparison, weather forecasts of temperature and precipitation are typically completely dominated by the seasonal trends after only a few days. Forecast error is lower for Oslo compared to Geilo and Tromsø , but relative to the average snow depth forecast error for Oslo is higher than for Geilo and Tromsø.
show forecasts for future snow depths in Oslo using model 2 for a season with little and much snow, respectively. The first, second and third row show forecasts for five, ten and three weeks into the future. In the left and the right column, we assume that zero and five days of reliable weather forecasts of temperature and precipitation are available, respectively. The solid curve shows the real snow depth data, while the dashed curves show 5% and 95% quantiles of the forecast distribution.
We see that given reliable weather forecasts more precise snow depth forecasts can be achieved. Comparing the bottom rows of Figures 8 and 9, 10 and 11 and 12 and 13, we see that the forecasts for three weeks into the future is quite different for seasons with little and much snow which shows that reliable long terms forecasts of snow depth are possible which are more precise than just using seasonal trends.
5 Closing remarks
This paper presents a first attempt to build statistical models for short and long term forecasts of snow depth. The results show that it is possible to do useful forecasts of snow depth long into the future. Further we found that model 2 (Section 3.2) perform better then model 1 (Section 3.2), but the advantage of model 1 compared to model 2 is that long term simultaneous scenarios of temperature, precipitation and snow depth is computed. This can be useful in for many applications. E.g. with respect to road safety the risk of slippery roads is especially high when the snow depth is above zero cm and at the same time the temperature is below zero C.
Several extensions to the suggested models are possible. Including other covariates like solar radiation, humidity, wind and the age of the snowpack may improve the forecasts (Roebber et al.; 2003; Kuusisto; 1984). Models that better separate sinking/aging from melting may be achieved by including the water content in the snow as a hidden layer in the model. Light snow tend to sink faster than denser snow and is not possible to separate in the model presented in this paper. Figures 1 3 reveal that the properties of snowfall and aging/melting varies at different locations. A natural extension is thus to include the models in this paper as part of a spatiotemporal model.
References
 (1)
 Blanchet and Davison (2011) Blanchet, J. and Davison, A. C. (2011). Spatial modeling of extreme snow depth, Annals of Applied Statistics 5(3): 1699–1725.
 Bocchiola et al. (2006) Bocchiola, D., Medagliani, M. and Rosso, R. (2006). Regional snow depth frequency curves for avalanche hazard mapping in central Italian Alps, Cold Regions Science and Technology 46(3): 204–221.
 Brown et al. (2003) Brown, R. D., Brasnett, B. and Robinson, D. (2003). Gridded North American monthly snow depth and snow water equivalent for GCM evaluation, AtmosphereOcean 41(1): 1–14.

Dunn and Smyth (2015)
Dunn, P. K. and Smyth, G. K. (2015).
dglm: Double Generalized Linear Models.
R package version 1.8.2.
http://CRAN.Rproject.org/package=dglm  Falarz (2002) Falarz, M. (2002). Longterm variability in reconstructed and observed snow cover over the last 100 winter seasons in Cracow and Zakopane (southern Poland), Climate Research 19(1): 247–256.
 Fnugg (n.d.) Fnugg (n.d.). fnugg.no. Accessed: 20161114.
 Jonas et al. (2009) Jonas, T., Marty, C. and Magnusson, J. (2009). Estimating the snow water equivalent from snow depth measurements in the Swiss Alps, Journal of Hydrology 378(1): 161–167.

Juga et al. (2013)
Juga, I., Nurmi, P. and Hippi, M. (2013).
Statistical modelling of wintertime road surface friction, Meteorological Applications 20(3): 318–329.
http://dx.doi.org/10.1002/met.1285 
Kienzle (2008)
Kienzle, S. W. (2008).
A new temperature based method to separate rain and snow, Hydrological Processes 22(26): 5067–5085.
http://dx.doi.org/10.1002/hyp.7131  Kohler et al. (2006) Kohler, J., Brandt, O., Johansson, M. and Callaghan, T. (2006). A longterm Arctic snow depth record from Abisko, northern Sweden, 1913 2004, Polar Research 25(2): 91–113.
 Kuusisto (1984) Kuusisto, E. (1984). Snow accumulation and snowmelt in Finland, Technical Report 55, National Board of Waters, Finland, Helsinki.
 McCabe and Wolock (2010) McCabe, G. J. and Wolock, D. M. (2010). Longterm variability in Northern Hemisphere snow cover and associations with warmer winters, Climatic Change 99(1): 141–153.
 Möller et al. (2013) Möller, A., Lenkoski, A. and Thorarinsdottir, T. L. (2013). Multivariate probabilistic forecasting using ensemble Bayesian model averaging and copulas, Quarterly Journal of the Royal Meteorological Society 139: 982–991.

R Core Team (2015)
R Core Team (2015).
R: A Language and Environment for Statistical Computing, R
Foundation for Statistical Computing, Vienna, Austria.
https://www.Rproject.org/  Roebber et al. (2003) Roebber, P. J., Bruening, S. L., Schultz, D. M. and Jr., J. V. C. (2003). Improving Snowfall Forecasting by Diagnosing Snow Density, Weather and Forecasting 18(1): 264–287.

Scherrer et al. (2013)
Scherrer, S. C., W thrich, C., CrociMaspoli, M., Weingartner, R. and Appenzeller, C. (2013).
Snow variability in the Swiss Alps 1864 2009, International
Journal of Climatology 33(15): 3162–3173.
http://dx.doi.org/10.1002/joc.3653  Sloughter et al. (2007) Sloughter, J. M., Raftery, A. E., Gneiting, T. and Fraley, C. (2007). Probabilistic quantitative precipitation forecasting using Bayesian model averaging, Monthly Weather Review 135(1): 3209–3220.
 Snow Forecast (n.d.) Snow Forecast (n.d.). snowforecast.com. Accessed: 20151109.

Stern and Coe (1984)
Stern, R. D. and Coe, R. (1984).
A Model Fitting Analysis of Daily Rainfall Data, Journal of
the Royal Statistical Society. Series A (General) 147(1): 1–34.
http://www.jstor.org/stable/2981736  Wen et al. (2013) Wen, L., Nagabhatla, N., Lu, S. and Wang, S.Y. (2013). Impact of Rain Snow Threshold Temperature on Snow Depth Simulation in Land Surface and Regional Atmospheric Models, Advances in Atmospheric Sciences 30(5): 1449–1460.
 Zhang et al. (2004) Zhang, Y., Li, T. and Wang, B. (2004). Decadal Change of the Spring Snow Depth over the Tibetan Plateau: The Associated Circulation and Influence on the East Asian Summer Monsoon, Journal of Climate 17(1): 2780–2793.
Comments
There are no comments yet.