Physical Modeling of Snow Gliding: A Case Study in the NW Italian Alps

: Snow gliding, a slow movement downhill of snow cover, is complex to forecast and model and yet is extremely important, because it drives snowpack dynamics in the pre-avalanching phase. Despite recent interest in this process and the development of some studies therein, this phenomenon is poorly understood and represents a major point of uncertainty for avalanche forecasting. This study presents a data-driven, physically based, time-dependent 1D model, Poli-Glide , able to predict the slow movement of snowpacks along a ﬂow line at the daily scale. The objective of the work was to create a useful snow gliding model, requiring few, relatively easily available input data, by (i) modeling snowpack evolution from measured precipitation and air temperature, (ii) evaluating the rate and extent of movement of the snowpack in the gliding phase, and (iii) assessing fracture (i.e., avalanching) timing. Such a model could be then used to provide hazard assessment in areas subject to gliding, thereby, and subsequent avalanching. To do so, some simplifying assumptions were introduced, namely that (i) negligible traction stress occurs within soil, (ii) water percolation into snow occurs at a ﬁxed rate, and (iii) the micro topography of soil is schematized according to a sinusoidal function in the absence of soil erosion. The proposed model was then applied to the “Torrent des Marais-Mont de La Saxe” site in Aosta Valley, monitored during the winters of 2010 and 2011, featuring different weather conditions. The results showed an acceptable capacity of the model to reproduce snowpack deformation patterns and the ﬁnal snowpack’s displacement. Correlation analysis based upon observed glide rates further conﬁrmed dependence against the chosen variables, thus witnessing the goodness of the model. The results could be a valuable starting point for future research aimed at including more complex parameterizations of the different processes that affect gliding.


Introduction
Snow gliding, defined as a gravity-driven, slow, and viscous downhill movement of snowpack, can lead to the formation of folds and cracks, which eventually may result in glide-snow avalanches [1]. The existence of glide processes has been recognized since the 1930s, and yet such a process is not fully understood given the lack of proper observations [2]. The forecasting of glide-snow avalanches contains major uncertainty, given their unpredictability [3][4][5][6]. Snow glide motion monitoring and modeling are essential for glide avalanche hazard assessment [4,7].
Gliding is very sensitive to the presence of free water at the snow-ground interface (e.g., [8,9]), and the formation of a soft, slushy soil film may largely influence the gliding 2 of 20 mechanism. Sources of free water at the snow-soil interface may include (i) rainfall; (ii) melt at the interface, resulting from heat storage in the soil; (iii) snowpack melt; and (iv) melt from geothermal hot spots and groundwater outflows (e.g., [10,11]), although the latter case is quite rare.
The current literature has focused on glide processes using a two-fold approach, namely (i) by investigating key parameters, such as water content, in the snow and at the snow-ground interface (e.g., [8,[11][12][13][14]); and (ii) by investigating the link between avalanche release activity and meteorological parameters (e.g., [4,6,8,15]). We also focused on the forces acting during gliding upon trees and structures [16].
Among others, Höller [1] recently provided a complete review of research conducted in the field of snow gliding and glide-snow avalanches, while also reviewing different approaches used since the 1930s, including descriptive methods, field measurements, models, and new technologies. Here, we propose a modeling approach, based upon flow measurements.
Among the first theoretical models, Heafeli [17] modeled the snow gliding process on an ideal plane, measuring the glide velocity against normal stresses on a glass plate at different inclinations. He found that the gliding velocity depends on the temperature and humidity of the glass plate (in the real world, the snow-ground interface). Most of the research has been based on a modeling approach, aiming to identify the relationship between the glide velocity and shear stress by considering the snow/ground driving features of gliding, such as soil temperature, surface roughness, slope and aspect, and the presence of free water at the interface (e.g., [18]). In der Gand and Zupanic [19] found that the glide velocity increases with the weight of snow and ground slope and decreases with increasing bottom friction. McClung [20] hypothesized and modeled three different mechanisms of snow gliding, namely (i) glide by creep, (ii) glide by regelation, and (iii) glide due to separation. Later, McClung and Clarke [21] proposed a model to take into account the role of free water at the snow-ground interface. Free water affects the glide velocity by reducing the snow's viscosity near the glide interface and by inducing partial separation of the snowpack from the interface. The formation of a wet basal layer favorable to snow gliding was studied, e.g., by Mitterer and Schweizer [11], who made a first attempt at modeling these processes. They underlined that the processes producing liquid water at the bottom of the snowpack can be divided, according to Clarke and McClung [8], into (i) cold temperature, and (ii) warm temperature events depending on whether the water at the base of the snowpack reaches from below (the soil) or from above (upper snow layers, [21]). Mitterer and Schweizer [11] created a 1D model to demonstrate that if dry snowpack is laid over a wet porous medium, the resulting large hydraulic pressure gradients may lead to an upward flux of the soil water content into the snowpack, a feature also found by Frigo et al. [22].
In addition to physical models of snow gliding, recent statistical models focused their attention on the most predisposing factors for intense snow gliding and glide-snow avalanche occurrence. Leitinger et al. [23] identified the key factors of forest stand, soil slope, winter precipitation, surface roughness, and aspect, and developed a spatial multiple regression model to produce maps of snow gliding in two areas in Austria. Other statistical models (e.g., [4,15]) have focused on the identification of driving factors for glide-snow avalanches by analyzing snow and weather in avalanche/no avalanche days. The main objective of this manuscript was to present a practical snow gliding model, Poli-Glide, requiring few, easily available input data. Specifically, we aimed to (i) reconstruct the snowpack evolution from measured precipitation and air temperature, with a special focus on the snow-soil interface; (ii) evaluate the rate and extent of the snowpack movement in the gliding phases; and (iii) assess the rupture (avalanching) time. We further carried out a correlation analysis to confirm the dependence of the observed glide rates against the chosen variables.

Study Area
The study area, located in Aosta Valley in the north-western Italian Alps, very close to the Mont Blanc Massif (4810 m asl), consists of an avalanche site called "Torrent des Marais-Mont de la Saxe", running on a west exposed slope, from 2115 m to 1250 m asl ( Figure 1). The long-term yearly mean precipitation in the study area is 840 mm yr −1 (1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010), and the mean annual air temperature is +2.8 • C (1993-2010, data from the AWS Courmayeur-Mont de la Saxe at 2076 m asl of the Ufficio Centro Funzionale, Aosta Region, [24]). The average cumulative annual snowfall is 275 cm at 1250 m asl  and about 450 cm at 2000 m asl [25]. The avalanche release area is typically characterized by intense snow gliding and by the formation of large glide cracks, often leading to the release of a glide-, often wet-, snow avalanche, mainly during springtime and sometimes in late autumn [18]. The crack/avalanche crown length usually ranges between 30 m and 100 m, depending on the intensity of the glide processes. On the left flank of the crack, a groundwater source is present. The area is characterized by a mean slope of 30 • , covered by pastures. Soil is frequently disturbed by the removal of the upper horizons and subsequent exposure of the subsoil, and signs of old and recent erosion, mainly due the snow-related processes, can be seen [12]. The bedrock is mainly made up of black clayish schists, calcareous sandstones, and, in some places, porphyritic granites.

Study Area
The study area, located in Aosta Valley in the north-western Italian Alps, very close to the Mont Blanc Massif (4810 m asl), consists of an avalanche site called "Torrent des Marais-Mont de la Saxe", running on a west exposed slope, from 2115 m to 1250 m asl ( Figure 1). The long-term yearly mean precipitation in the study area is 840 mm yr −1 (1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010), and the mean annual air temperature is +2. 8 °C (1993-2010, data from the AWS Courmayeur-Mont de la Saxe at 2076 m asl of the Ufficio Centro Funzionale, Aosta Region, [24]). The average cumulative annual snowfall is 275 cm at 1250 m asl  and about 450 cm at 2000 m asl [25]. The avalanche release area is typically characterized by intense snow gliding and by the formation of large glide cracks, often leading to the release of a glide-, often wet-, snow avalanche, mainly during springtime and sometimes in late autumn [18]. The crack/avalanche crown length usually ranges between 30 m and 100 m, depending on the intensity of the glide processes. On the left flank of the crack, a groundwater source is present. The area is characterized by a mean slope of 30°, covered by pastures. Soil is frequently disturbed by the removal of the upper horizons and subsequent exposure of the subsoil, and signs of old and recent erosion, mainly due the snow-related processes, can be seen [12]. The bedrock is mainly made up of black clayish schists, calcareous sandstones, and, in some places, porphyritic granites.

Data
Snow and weather data were provided by the automated weather station Pré-Saint-Didier Plan Praz (2044 m asl, 8.5 km away from the avalanche release area but with a similar altitude and the same aspect), managed by the Ufficio Centro Funzionale (UCF)-Aosta Region. Data from the Courmayeur Mont de La Sax station, closer to the study site, could not be used. Therefore, snow depth data were not dependable, since large-scale erosion, or (wind) accumulation phenomena occur. The station of Pré-Saint-Didier Plan Praz has been active since 2002, and it considered well representative of the area, especially for the snow depth at 2000 m asl. The parameters available from this AWS were air temperature ( • C), snow depth (cm), rain (mm), and solar radiation (Wm −2 ). To assess the physical properties of the snow cover, several snow profiles, taken according to the procedure in Fierz et al. [26], were dug close to the release area in a safe place, where avalanches rarely occur.
In the release area of the glide-snow avalanches, specific instrumentation was installed in summer 2009 to continuously measure snow gliding and some physical properties of snow and soil. A detailed description of the experimental set up is given in the work of Ceaglio et al. [12]. Here, we report the main information regarding the collected data, namely snow gliding (cm), using snow glide shoes placed at the snow/soil interface, snow temperature and water content in the basal layer of the snowpack, and Soil temperature and water content at 0 cm (snow/soil interface), −5 cm, −15 cm into the soil. Data loggers were set up to scan measures of different parameters every minute and to store average values every 30 min.

Poli-Glide Model
Here, we developed a physically based snow gliding (SG) model, Poli-Glide. The proposed model is mono-dimensional (1D), distributed (as per a number of cells along a main flow line), physically based, and works at a daily level. This model builds upon a number of modules, each one simulating a specific snow-related process among those more relevant for the depiction of gliding. Gliding modeling requires the depiction of a number of processes, starting from snowfall and accumulation, to subsequent snowpack dynamics, energy and water balance, and to subsequent creeping motion, or gliding. Here, these processes are schematized as (i) snow accumulation, (ii) snowpack settling, (iii) snowpack ablation, (iv) snowpack temperature distribution, (v) snowpack and soil water content and budget, and (vi) snowpack gliding. Each of these phenomena, or processes cascading into one another, is modeled according to the available approaches within the present literature, and their combined use leads to the final Poli-Glide model. We depict our approach to each process below.

Snow Accumulation
Each time precipitation, P, occurs, partial or total snowfall, P s (and complementary liquid precipitation, P l ), is predicted against air temperature, T a , as with α l , the fraction of liquid water, depending upon temperature and T inf and T sup are defined/tuned against observations. Here, we used T inf = −0.5 • C and T sup = 0.5 • C. When precipitation occurs and T a ≤ T sup , solid precipitation (snowfall) with T a expressed in • C (e.g., for T a = 0 • C, ρ 0 = 152.5 kgm −3 ). Here, a density lower than 50 kgm −3 cannot be calculated. Few fresh snow measurements were carried out in the surroundings of the glide shoes, and they all displayed values higher than 50 kgm −3 . Additionally, as reported by Bocchiola and Rosso [28], who reported 2500+ measurements of fresh snow within the central Italian Alps, fresh snow density was very rarely (8 cases) below 50 kgm −3 . Accordingly, Equation (3) would be suitable most of the time in practice.

Snowpack Settling
As snowpacks evolve in time during winter, snow compaction occurs, with the decrease in snow depth, H s (e.g., [29]). Settling in time for one single snow layer may be modeled, according to a power law, as with H 0 being the initial depth of the new snow, n the number of days after a snowfall event (day 0), and m an exponent to be tuned against snow data (e.g., m ≈ 0.3, [29]). Neglecting mass losses (i.e., for negligible snow melt, and sublimation) and snow redistribution by wind, one may estimate settled snow density by assuming mass conservation as Given the superposition of different layers as given by subsequent snowfalls, settling of each snow layer may accelerate. However, here, daily changes in snow density are tracked by the model for each and every snow layer (and one layer i for each of the N s events until day j) independently. At the end of each day j, the mean snowpack density, ρ SP,j , of the snowpack, H S,j , is taken as

Snowpack Ablation
Here, snow ablation is simulated according to a degree-day approach, able to sufficiently capture the process at a daily scale [30,31]. In doing so, snow melt, M s , is as follows: with being a melting factor to be tuned against data, and T t a threshold temperature, often set to T t = 0 • C, as we did here after a preliminary screening of snowpack data.

Snowpack Temperature Distribution
Temperature within the snowpack is assumed to vary linearly (i.e., along a downward coordinate z starting from snow surface), as a function of the temperature of snow surface at the interface with air (T s,a ), and the temperature of snow at the bottom at the interface with soil (T s,s ), as follows: with temperature changing along the snow depth as Here, snow to soil interface temperature, T s,s , was measured by way of temperature probes, as described. In our model, T s,s was used to assess potential gliding (i.e., when T s,s = 0 • C). Whenever T s,s was not available, one could use T a to assess the temperature profile in the snowpack (and T s,s ), say, by way of calculation of the energy budget within the snowpack as driven by weather patterns (i.e., solar radiation, wind, air moisture, etc.).

Snowpack and Soil Water Content and Budget
Water storage at the snow-soil interface is modeled according to the water budget of soil, or soil water content, S (expressed here in mm, i.e., as a water volume over a unit surface area), as conditioned upon snow melt (e.g., [32]), as follows: with S j [mm] soil water storage at day j, depending upon the content the day before, S j−1 , liquid precipitation in the same day, P l,j , snow melt in the same day, M s,j , and soil infiltration in the same day, Q g,j (the three latter variables in mmd −1 , and ∆t time step, 1 day here) [33]. For simplicity, here, is assumed that snow melt water percolates within the snowpack and reaches the ground within one day. Underground infiltration is as follows [33]: with K sat [mmd −1 ] being the saturated soil conductivity and S max [mm] the largest (potential) soil water storage. Aboveground water depth at ponding (i.e., in the case of saturated soil the depth of water flow reached over the surface), h w , is then

Gliding Velocity
Here, we rely upon a gliding velocity formulation by McClung and Clarke [21], linking gliding velocity, U g [ms −1 ], to shear stress, τ [Pa], at the snow-soil interface With D* [m] being the depth of stagnation depth, µ [Pas] the shear viscosity, and ν [.] the Poisson snow viscosity coefficient. Figure 1 reports a sketch of the line where the approach by McClung and Clarke [21] is applied. The key variable here is D*, depending upon soil geometry, and water depth at the snow-soil interface. If one models soil surface roughness using a sinusoidal function with wavelength, λ 0 , and amplitude, A [16,20,21,34], in the presence of a water layer, h w , stagnation depth, D*, can be estimated as Whenever h w = A, soil roughness is null, and D* would no longer be defined. Equation (14) represents an extension of the one provided by Salm [34], based upon further elaboration by McClung and Clarke [21]. Therein, a partial drowning of the bed, with roughness defined using a sinusoidal function, contributes to a decrease in soil drag capac-ity, and thereby to an increase in the stagnation depth, D * . Shear stress in the presence of a snowpack, H S , with density, ρ S , is taken from the balance of gravity and friction forces along the slope as with φ [ • ] being the snow-soil friction angle and c S [Pa] the snowpack cohesion. Accordingly, gliding velocity may be taken as depending upon snowpack depth, density, mechanic properties (cohesion and friction angle), and soil roughness.

Model Set Up
The gliding model was set up for the case study area under a 1D distributed scheme ( Figure 1). Gliding motion occurs along a pre-defined flow line, according to our experimental set up. We sketched four different flow lines, two for the A plot (A1 and A2 glide shoes) and two for the B plot (B1 and B2 glide shoes).
Using a shape file depicting the 2010 avalanche event and uploading the coordinates of the glide shoes within the GIS tool, we found out the release zone of the potential gliding movements. Then, using the iso-level curves of the slope, we depicted the flow lines in the snow flowing area. Each flow line was then discretized into a number of cells with 2 × 2 m 2 sides. For each cell in each flow line, we then estimated altitude and slope. Upstream cells of each line are at ca. 2118 m a.s.l., and downstream cells are at ca. 2035 m a.s.l., with each line having an average slope of ca. 29 • . For the purpose of setting up the gliding model for the Mont de la Saxe site, a number of steps were carried out, as follows.
Precipitation data were gathered and elaborated to obtain snowfall (i.e., snow accumulation) estimation. The AWS stations of Prè-Saint-Didier Plan Praz and Courmayeur Mont de la Saxe deployed unheated pluviometers, unable to gather new snow during snowfall events. Accordingly, during winter, very low amounts of precipitation were recorded. To properly assess the fresh snow amount, we proceeded as follows. Using the snowpack depth, H S , data from Prè-Saint-Didier Plan Praz (winter 2010-2011), we deduced snowfall events ex post from the observed peaks of snowfall depth, or ∆H S > 0. Given the known air temperature, T a , we assessed snow density via Equation (3) and subsequently estimated the new snow mass (i.e., new snow water equivalent SWE 0 ) as DH S i.e., trading new snow depth for the snowpack depth variation (in day 0). We subsequently estimated the Martinec exponent, m, in Equation (4). In the absence of snow melt (i.e., when Equation (5) is valid), only snow settling occurs, and m could be estimated by iterations, until snowpack settling (i.e., decrease in snowpack depth H S during dry days) at Prè-Saint-Didier Plan Praz could be well mimicked (winter 2010-2011) using Equation (4). The melt factor in Equation (7) was also estimated iteratively from snow data, by modeling snowpack depth during melting events (i.e., days with T a > T t = 0) until full snowpack depletion.
Thermal lapse rates during the snow season were assessed on a monthly basis, and independently for year 2010 and 2011, using the data from 19 temperature stations (including Prè-Saint-Didier Plan Praz and Courmayeur Mont de la Saxe) belonging to the UCF of Val d'Aosta, ranging from 935 m a.s.l. to 2430 m a.s.l. Notice that the considered gliding sites span over a relatively narrow range of altitude (2035 m a.s.l. to 2118 m a.s.l. as reported); thus, temperature changes would not be so large in a first approximation. However, given the effect of temperature upon snow melt and water accumulation within the snowpack Climate 2021, 9, 171 8 of 20 and at the soil-snow boundary, we tried to accurately depict thermal changes with altitude. The obtained lapse rate varied from −2.1 • Ckm −1 during January 2010 to −6.7 • Ckm −1 during April 2011.
The hydraulic conductivity of soil below snowpacks depends upon soil composition. Here, soil texture was found to mainly feature sand (81.6% in weight) with some silt (16.2%) and very little clay (2.2%), thus pointing to a rather permeable soil, with an estimation of K sat = 234 mmd −1 [18]. Considering an average soil depth of ca. 400 mm as from local soil maps and a saturation coefficient of ca. 31% as from soil texture [18], one achieved an estimated maximum soil storage of S max = 124 mm. The values of shear viscosity and the Poisson snow viscosity used here were µ∼5 × 10 10 Pas, and ν∼0.2, as reported by Ancey and Bain [35], for typical values of this variable for snow. Gliding velocity parameters (λ 0 , A, φ, and c S , Equations (15) and (16)) are difficult to be obtained in practice, and here, we used them as tuning parameters to be assessed against gliding observations.

Correlation Analysis
We carried out a correlation analysis, considering the daily values of glide rate, i.e., displacement of each day as a dependent (predictable) variable. Some potentially predictive (independent) parameters were taken, e.g., snow cover, depth and temperature, air temperature, precipitation, solid and liquid, and the interface temperature and water content in the soil, which were relevant, as reported in literature (e.g., [13]).

Snowpack Dynamics
Albeit our model was not developed explicitly to simulate snowpack dynamics, and more accurate models exist for the purpose [36]. The state of snowpacks influences gliding, and therefore, reasonably accurate mimicry of the former is necessary for proper modeling of the latter. Therefore, here, we analyze how snowpack dynamics are depicted by our model after parameter tuning.
First, we estimated the snow melting factor using Equation (7) Table 1.
The tuning (calibration) value was m = 0.13 with acceptable statistics (R 2 = 0.96, see Bias % and RMSE in Table 1). When validation was carried out, some loss of accuracy was seen (R 2 = 0.52), but the results seemed still acceptable. The obtained value of Martinec's exponent could therefore be taken as a constant during the winter seasons, without a large loss of accuracy. The modeled values of H S quite accurately mimic the observed ones in the accumulation phase, where new snowfall and settling are the prevailing phenomena. Our approach is thus able to accurately assess snowfall occurrence (i.e., snowfall events), snow accumulation and subsequent compaction of the snow crust, including the combined effect of different snow layers.

Snow Gliding Parameterization
The four parameters (cohesion factor, angle of friction, wavelength, and the amplitude of terrain topography) that best mimicked the measured displacements were then identified. The calibration of these parameters was necessary here because experimental measurements were not available, however difficult in the field, and given the lack of literature information, this made the assessment uncertain. To make up for such uncertainty, several simulations were then carried out. We firstly varied the cohesion factor and the friction angle, the same for each glide shoe, and subsequently the wavelength and amplitude of ground topography, tuning the values for each glide shoe separately. Regarding the topography of terrain, Höller et al. [16] proposed that for a generic mountain slope, the wavelength would generally be less than 2 m or so, while amplitude may be between 10 and 500 mm. However, no correlations were proposed (depending, for example, on the type/texture of soil) which could be useful to estimate such parameters. Otherwise, in their application, Ancey and Bain [35] found a wavelength of 10 m, outside the range found by Höller et al. [16], and an amplitude of 100 mm. However, these parameters are rather site specific. The local assessment of these values is needed, but this requires specific field procedures that are not currently available for the Mont de La Saxe site.
Regarding cohesion and angle of friction, In der Gand and Zupancic [19] found a value of friction of 0.3 (i.e., a friction angle φ = 16.7 • ). Podolskiy et al. [37] showed that the best combination for the evaluation of snow stress according to a Mohr Coulomb criterion is a cohesion of 1.6 kPa, with a variable friction angle between φ = 22.5 • and 60 • , a rather large range. Thereby, a wide range of possible values for φ was explored in tuning our model. However, after several simulations, we observed that a constant angle of friction would not provide accurate results at all sites. Instead, by considering an angle of friction changing against cumulative displacement of the snowpack, we found a much better depiction of the gliding process. In particular, we found a decrease in the friction angle as a function of displacement. Such a circumstance may be due to changes in roughness of the ground surface, by flattening of soil. Additionally, since the permeability of the soil can be much lower than that of snow, a layer of water, resulting from melting (and excess of soil moisture below), can persist over the soil-snow interface, reducing the friction of soil, as found, e.g., by [12]. In that study, the authors identified the liquefaction of soil as a potential factor capable of contributing to the processes of snow gliding, as well as soil erosion, through the formation of a layer of mud and water at the soil-snow interface, capable of reducing roughness.
It was therefore decided to include an angle of friction changing within a given range in the model, i.e., decreasing linearly from an initial maximum value, down to a minimum one, which then would remain constant until the end of simulation. This approach, compared to a constant or sudden (ramp-like) change of the angle friction, brought better results in the tuning phase, and it was therefore chosen to describe gliding. In Figure 3, we present the results for the four glide shoes for season 2010. We report therein the modeled and measured snow gliding, the measured temperature of air and soil, measured precipitation, and snow depth modeled over time, as described above.  Glide measurements in Figure 3 show that the snowpack slipped differently in two glide shoes, A1 and A2, albeit they were in close range. The gliding of snowpacks may indeed be strongly influenced by local topography, and even small changes can cause different movements. Tracing of the flow lines was carried out, starting from a digital model of the terrain, with a resolution of 2 m, and with the help of slope contours. Since it was not possible to accurately reconstruct the paths followed by the two sleds, and since these were very close together, here, we also decided to mimic the movement of a virtual glide shoe, placed between A1 and A2, which we call glide shoe A12, and displaying a movement composed as a mean of the movements of A1 and A2 glide shoes. In Table 2, we report the values of the parameters chosen for best modeling of the phenomenon. For each glide shoe, a value of the wavelength, λ0, and amplitude, A, of the terrain topography was found. In spite of some observed variation among glide shoes, somewhat consistent values seem to occur. Moreover, the values are comparable to those found by Ancey and Bain [35] and are larger than those found by Höller [1] and Leitinger et al. [23], albeit within the same order of magnitude.

Gliding Events
The remaining parameters were taken as equal for all glide shoes, as seen in Table 2, given that field surveys [12] displayed a substantially constant soil composition in the area of interest. Glide measurements in Figure 3 show that the snowpack slipped differently in two glide shoes, A1 and A2, albeit they were in close range. The gliding of snowpacks may indeed be strongly influenced by local topography, and even small changes can cause different movements. Tracing of the flow lines was carried out, starting from a digital model of the terrain, with a resolution of 2 m, and with the help of slope contours. Since it was not possible to accurately reconstruct the paths followed by the two sleds, and since these were very close together, here, we also decided to mimic the movement of a virtual glide shoe, placed between A1 and A2, which we call glide shoe A12, and displaying a movement composed as a mean of the movements of A1 and A2 glide shoes. In Table 2, we report the values of the parameters chosen for best modeling of the phenomenon. For each glide shoe, a value of the wavelength, λ 0 , and amplitude, A, of the terrain topography was found. In spite of some observed variation among glide shoes, somewhat consistent values seem to occur. Moreover, the values are comparable to those found by Ancey and Bain [35] and are larger than those found by Höller [1] and Leitinger et al. [23], albeit within the same order of magnitude.
The remaining parameters were taken as equal for all glide shoes, as seen in Table 2, given that field surveys [12] displayed a substantially constant soil composition in the area of interest. Table 2. Model parameters of the best-fitting model. The topography of the terrain values or each of the 4 glide shoes (A1, A2, B1, B2) and the virtual one (A12). The cohesion and the angle of friction, as parametrized in paragraph 5.2, are the same for each glide shoe. We evaluated the objective goodness of fit using Nash-Sutcliffe NSE efficiency Root Mean Square Error, RMSE (e.g., [38]), of cumulative glide (in cm), and glide rate (in cm/d).

A1
In Table 3, we also report the maximum (absolute) difference between the modeled and measured values of cumulative glide. Figure 3 shows the trend of the measured and modeled cumulative movement of snow. Table 3. Goodness of fit statistics for snow glide, rate, and cumulative glide for each of the 4 glide shoes (A1, A2, B1, B2) and the virtual one (A12). Model efficiency after calibration is slightly better for the B glide shoes than for the A ones (and especially low for glide shoe A1, see Figure 3a).

A1
However, the gliding model seems to provide acceptable results, indicating that snow gliding can be modeled acceptably, at least in a first approximation, with the Poli-Glide model. For the glide shoes A1 and A2, during winter 2010, there is a good correspondence between the measured and modeled final cumulated glide, but less correspondence exists between glide rate in the central period. If one considers the virtual glide shoes, A12 (Figure 3e), a much better performance is obtained with a greater correspondence between the modeled and (average) observed displacement of the two glide shoes. This circumstance may be indicative of a dependence of the glide phenomenon upon local (micro)topographic characterization, somewhat filtered out by averaging between the two glide shoes, or possibly of the effect of some noise in measurements. Instead, the patterns modeled for glide shoes B1 and B2 show a substantial correspondence with the measured values.

Winter 2011
The tuning parameters (λ 0 , A, φ, and c S ) that best simulate the trend of displacements measured during the 2010 season were initially also used for the simulations of gliding during the 2011 winter season (reported in [12], Figure 9b). For the 2011 season, data are available from late autumn (8 November 2010) to late spring (30 April 2011). No movement was recorded from the A1 glide shoe ( [12], Figure 9b). Albeit the model could capture snow mantle dynamics well during winter 2011 (Figure 2b), and tracking gliding using the same values of the parameters as tuned for the 2010 season led to less satisfactory results ( Figure 4). Using the same values of cohesion and friction angle, the model still managed to identify the (timing of) movement around mid-January, but it later deviated and failed to describe the final cumulative displacement recorded. The date when the model started deviating more corresponds to 17 January 2011, when a daily movement of about 1.4 m was recorded for all glide shoes (A2, B1, B2), i.e., more than 50% of the total movement, which indeed may not be considered as a typical snow gliding phenomenon (see [12], Figure 9b).
The unsatisfactory results obtained by applying the model to the 2011 winter season (after the calibration of the parameters on the 2010 season) can be likely charged to the substantial differences found between the two winter seasons, from the point of view of both displacement of snowpack dynamics and the weather.
The movements recorded by the glide shoes during the two winter seasons in 2010 and 2011 are very different from each other both concerning the trend over time and the cumulative movement at the end of the season. These two hydrological years were indeed characterized by different meteorological conditions [12]. During the 2011 winter season, there were no significant correlations between the movement of the snowpack and the various meteorological, snow, and pedological variables. However, the few but intense snow slips were qualitatively associated with significant settling of the snowpack and rises in air temperatures [12]. Further, the presented Poli-Glide model was cast so as to (satisfactorily) simulate the gliding phenomenon as long as it occurs as a continuous (quasistatic) movement. Difficulties are instead encountered in days when large movements are recorded, associated with the possible opening of cracks. Snow gliding and the formation of fractures are indeed different phenomena, characterized by different distributions of stresses. The snow gliding formulations introduced in the model indeed consider slower snow movements and generally lower daily glide rates, and the model was therefore unable to simulate more impulsive events as tracked in winter 2011. capture snow mantle dynamics well during winter 2011 (Figure 2b), and tracking gliding using the same values of the parameters as tuned for the 2010 season led to less satisfactory results (Figure 4). Using the same values of cohesion and friction angle, the model still managed to identify the (timing of) movement around mid-January, but it later deviated and failed to describe the final cumulative displacement recorded. The date when the model started deviating more corresponds to 17 January 2011, when a daily movement of about 1.4 m was recorded for all glide shoes (A2, B1, B2), i.e., more than 50% of the total movement, which indeed may not be considered as a typical snow gliding phenomenon (see [12], Figure 9b).
The unsatisfactory results obtained by applying the model to the 2011 winter season (after the calibration of the parameters on the 2010 season) can be likely charged to the substantial differences found between the two winter seasons, from the point of view of both displacement of snowpack dynamics and the weather.
The movements recorded by the glide shoes during the two winter seasons in 2010 and 2011 are very different from each other both concerning the trend over time and the cumulative movement at the end of the season. These two hydrological years were indeed characterized by different meteorological conditions [12]. During the 2011 winter season, there were no significant correlations between the movement of the snowpack and the various meteorological, snow, and pedological variables. However, the few but intense snow slips were qualitatively associated with significant settling of the snowpack and rises in air temperatures [12]. Further, the presented Poli-Glide model was cast so as to (satisfactorily) simulate the gliding phenomenon as long as it occurs as a continuous (quasi-static) movement. Difficulties are instead encountered in days when large movements are recorded, associated with the possible opening of cracks. Snow gliding and the formation of fractures are indeed different phenomena, characterized by different distributions of stresses. The snow gliding formulations introduced in the model indeed consider slower snow movements and generally lower daily glide rates, and the model was therefore unable to simulate more impulsive events as tracked in winter 2011.
(a)    Figure 5 reports (in a polar chart format) the results of the correlation analysis for each of the four gliding shoes. A strong positive correlation of gliding speed, U, can be seen against water content in the soil and against snow depth. Indeed, these are the main parameters used by the Poli-Glide model to assess gliding speed. Within the model, we assumed the dependence of glide speed on the depth of snowpack (H S ), and even upon stagnation depth (D * ), with the latter being a function of the water content in the soil in practice. These results agree with those obtained, e.g., by Ceaglio et al. [13], where the authors analyzed correlations of gliding rate against a number of variables in detail, evaluating use of empirical gliding models. Accordingly, the proposed Poli-Glide model seemingly considers proper drivers for the depiction of snow gliding. Figure 5 reports (in a polar chart format) the results of the correlation analysis for each of the four gliding shoes. A strong positive correlation of gliding speed, , can be seen against water content in the soil and against snow depth. Indeed, these are the main parameters used by the Poli-Glide model to assess gliding speed. Within the model, we assumed the dependence of glide speed on the depth of snowpack ( ), and even upon stagnation depth ( * ), with the latter being a function of the water content in the soil in practice. These results agree with those obtained, e.g., by Ceaglio et al. [13], where the authors analyzed correlations of gliding rate against a number of variables in detail, evaluating use of empirical gliding models. Accordingly, the proposed Poli-Glide model seemingly considers proper drivers for the depiction of snow gliding. Figure 5. Correlation analysis of gliding speed, , against chosen variables: air temperature, , interface temperature, , , soil water content, , solid precipitation, , liquid precipitation, , total precipitation, , snow depth, and snow temperature, .

Conclusions
Snow gliding is very difficult to predict and model. In this study, we tried to contribute to the understanding of this phenomenon by presenting a physically based and time-dependent one-dimensional snow gliding model, able to describe the movement of a snowpack along a flow line on a daily scale.
We were able to create a gliding model that required the least amount of input data. The model can sufficiently reconstruct the evolution of a snowpack, at least under a viscous flow condition. Starting from measured precipitation values, the model can track snow thickness, density, thermal conductivity, SWE, and temperature for each of the individual layers in the snow mantle. The model can also track snow velocity and displacements caused due to snow gliding. Therefore, the model provides a potentially useful tool in assessing the risks in areas subject to these phenomena or in estimating the impacts therein, such as soil erosion. To allow the modeling, simplifying hypotheses were introduced. Namely, daily assessment was carried out, which may have neglected processes occurring at a shorter time scale (e.g., hourly). Water percolation through a snowpack resulting from precipitation or melting was not considered, and it was assumed that percolating water reaches the ground-snow interface within a day. Figure 5. Correlation analysis of gliding speed, U, against chosen variables: air temperature, T a , interface temperature, T s,s , soil water content, S, solid precipitation, P s , liquid precipitation, P l , total precipitation, P, snow depth, H S and snow temperature, T S .

Conclusions
Snow gliding is very difficult to predict and model. In this study, we tried to contribute to the understanding of this phenomenon by presenting a physically based and timedependent one-dimensional snow gliding model, able to describe the movement of a snowpack along a flow line on a daily scale.
We were able to create a gliding model that required the least amount of input data. The model can sufficiently reconstruct the evolution of a snowpack, at least under a viscous flow condition. Starting from measured precipitation values, the model can track snow thickness, density, thermal conductivity, SWE, and temperature for each of the individual layers in the snow mantle. The model can also track snow velocity and displacements caused due to snow gliding. Therefore, the model provides a potentially useful tool in assessing the risks in areas subject to these phenomena or in estimating the impacts therein, such as soil erosion. To allow the modeling, simplifying hypotheses were introduced. Namely, daily assessment was carried out, which may have neglected processes occurring at a shorter time scale (e.g., hourly). Water percolation through a snowpack resulting from precipitation or melting was not considered, and it was assumed that percolating water reaches the ground-snow interface within a day.
The results also show that further investigations of the phenomenon are needed to better determine the contribution of the various factors and parameters used in the calibration of the model. In particular, it would be necessary to deepen and improve the definition of the wavelength and the amplitude of the terrain topography, investigating their size more accurately through field surveys, and possibly linking their values to soil textures and cover (e.g., vegetation). Further studies must be carried out to link cohesion and friction to soil characteristics for practical determination.
One other possible improvement could be in the spatial analysis of the phenomenon with the use of a remote sensing tool able to evaluate snow cover areas [39], or depths [40] but also able to evaluate the snow water equivalent [41]. In this way, one can obtain the spatial distribution characteristics of the snow mantle and is able to apply the proposed methodology (or a similar one) independently of the presence of sensors.
More rapid snow movements, as tracked during winter 2011, and possibly due to cracking of snow mantle, could not be modeled using our approach, and further effort is required in this direction.
Our Poli-Glide model can be seen as a good starting point for future investigations aimed at expanding knowledge of the processes affecting slow gliding, and furthermore, this model can be improved to be applied in case of ice avalanches or rocky ice avalanches.
The expected changes in the behavior of the cryosphere under impending climate change hereon include potential for the more frequent occurrence of wet snow precipitation and rain on snow events, even at high altitudes, resulting in a general evolution of the snowpack towards warmer and wetter conditions [4]. As such, the occurrence of meteorological conditions favorable to triggering of avalanches due to gliding snow and/or flowing/cracking ice is likely to increase, as visible from a recent calamitous event [42]. The monitoring and modeling of snow gliding, and possibly even of crack opening in snow and ice, including with use of Poli-Glide, could represent an essential step towards managing snow/ice avalanche hazards and associated risks [4].