A local model of snow–firn dynamics and application to the Colle Gnifetti site

The regulating role of glaciers in catchment run-off is of fundamental importance in sustaining people living in low-lying areas. The reduction in glacierized areas under the effect of climate change disrupts the distribution and amount of run-off, threatening water supply, agriculture and hydropower. The prediction of these changes requires models that integrate hydrological, nivological and glaciological processes. In this work we propose a local model that combines the nivological and glaciological scales. The model describes the formation and evolution of the snowpack and the firn below it, under the influence of temperature, wind speed and precipitation. The model has been implemented in two versions: (1) a multi-layer one that considers separately each firn layer and (2) a single-layer one that models firn and underlying glacier ice as a single layer. The model was applied at the site of Colle Gnifetti (Monte Rosa massif, 4400–4550 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$). We obtained an average reduction in annual snow accumulation due to wind erosion of 2×10³ $\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ to be compared with a mean annual precipitation of about 2.7×10³ $\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$. The conserved accumulation is made up mainly of snow deposited between April and September, when temperatures above the melting point are also observed. End-of-year snow density, instead, increased an average of 65 kg m⁻³ when the contribution of wind to snow compaction was added. Observations show a high spatial and interannual variability in the characteristics of snow and firn at the site and a correlation of net balance with radiation and the number of melt layers. The computation of snowmelt in the model as a sole function of air temperature may therefore be one of the reasons for the observed mismatch between model and observations.

1 Introduction

Glacier ice covers almost 16×10⁶ km² of the Earth's surface, of which it is estimated that only 3 % is retained by the mountains outside the polar regions (Benn and Evans, 2010). Despite this small percentage, the amount of water stored in mountain glaciers plays a key role in sustaining people living in low-lying areas (Adhikary, 1993), influencing run-off on a wide range of temporal and spatial scales (Jansson et al., 2003; Huss et al., 2010). Storing water coming from precipitation in winter and delaying the time in which it reaches the river network, mountain glaciers sustain streamflow in hotter and drier periods when precipitation is lacking and when it is most needed for agriculture and as drinking water (Fountain and Tangborn, 1985; Hagg et al., 2007).

Jost et al. (2012) studied a Canadian river basin, covered for only 5 % by glaciers, and they found that ice melt contributes up to 25 % to streamflow in August, up to 35 % to streamflow in September, and between 3 % and 9 % to total streamflow.

In high mountain river basins of the northern Tian Shan (central Asia), with areas of glaciation higher than 30 %–40 %, glacier melt contribution is 18 %–28 % of annual run-off but it can increase to 40 %–70 % during summer (Aizen et al., 1996).

The reduction in glacier volume observed over the past 150 years (Vaughan et al., 2013; Hock et al., 2019) will result in a change in the present distribution and amount of water storage and release, with implications for all aspects of watershed management (Hock et al., 2005) and with consequently high economic impacts (Huss et al., 2010). The prediction of these changes is therefore fundamental in order to assess and reduce their impacts, optimizing consequently the management of water resources. To accomplish this task, models that integrate hydrological, nivological and glaciological components and that consider a variable glacier extension and the transient response of glaciers to climate change are required (Luo et al., 2013).

Despite their importance, fully integrated glacio-hydrological catchment models are not common in the literature (Wortmann et al., 2019). Some examples of glacio-hydrological models are provided by the works of Huss et al. (2010), Naz et al. (2014), Seibert et al. (2018) and Wortmann et al. (2019).

Wortmann et al. (2019) grouped the main problems of glacio-hydrological models into two categories: integration and scale. With integration problems they refer to the simplified or absent description of the remaining catchment hydrology in models that describe glacier processes in detail. The decrease in the fraction of ice-covered areas requires a proper description of both components, even in basins that are currently highly glacierized. Another aspect is the integration of nivological and glaciological components: a joint simulation of glacier mass balance and snow accumulation and melt is required in order to avoid inconsistencies (Jost et al., 2012; Naz et al., 2014). The problems of scale arise from the different resolutions required by glacial, nivological and hydrological processes. Physically based models that consider all glacier processes (mass balance, subglacial drainage and ice flow dynamics) are often too computationally expensive to be used in a combined glacio-hydrological model that considers the entire catchment. In addition, they are characterized by a complexity higher than the one of many semi-distributed hydrological models. It is therefore necessary to develop glacier models with a degree of complexity similar to the one of hydrological models but that are still able to reproduce important processes (Seibert et al., 2018).

In the present work, we give our contribution proposing a local model that follows the transformation of snow into firn and glacier ice under the influence of meteorological variables (temperature, precipitation and wind speed). Existing firn densification models are, in general, forced by snow characteristics. In this sense, the presented model allows us to move the boundary of the firn densification models from surface accumulation and density to hourly meteorological series. When we do not assume stationarity in the climate, in fact, this is required to properly capture the effects of climate changes.

The core of the model was derived from mass balance, momentum balance and rheological equations, governing the evolution of snowpack and firn (depth and density of snow and firn, depth of water and refrozen meltwater and rain inside the snowpack). The calculations of the terms in the resulting equations were then approached looking at methods already used in the literature; for example, snowmelt mass flux was computed with a temperature-index approach and the run-off from the snowpack with a flow matrix approach. We present two versions of the model: (1) a version (multi-layer) that considers separately each firn layer and (2) a version (single-layer) that models firn and underlying glacier ice as a single layer. The latter consists of only six equations, and it is therefore more suitable for a possible application to a hydrological model. The former consists of four equations for the snowpack plus two equations for each firn layer. Providing a profile of density with depth, it captures better the influence of meteorological variables on snow and firn characteristics. Besides, it allows a better validation of the snow–firn model. The equations that describe the snowpack are derived from the work of De Michele et al. (2013) and later Avanzi et al. (2015), modified in order to take into account the contribution of wind erosion and the transformation of snow into firn. To model the firn component, both the densification model of Arnaud et al. (2000) and the one of Herron and Langway (1980) were implemented. In order to test the model, a high-altitude site, Colle Gnifetti, belonging to the Monte Rosa massif, was chosen.

The paper is organized as follows: we present the model in Sect. 2, illustrate the case study in Sect. 3, give the results in Sect. 4 and discuss them in Sect. 5. The conclusions are given in Sect. 6.

2 Methodology

In this section, firstly the snowpack model, proposed by De Michele et al. (2013) and later modified by Avanzi et al. (2015) with the addition of the contribution of wind to snow transport, is illustrated and secondly the model with the integration of snow and firn processes is presented.

2.1 Snow model

The snowpack is modelled, according to De Michele et al. (2013) and Avanzi et al. (2015), as a mixture of dry and wet constituents. The solid deformable skeleton that consists of both snow grains and pores has a total volume V_S, unit area, height h_S, mass M_S and density ρ_S. The liquid water inside the pores has a volume V_W, unit area, height h_W, mass M_W and constant density ${\mathit{\rho }}_{\mathrm{W}}=\mathrm{1000}\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$. The refrozen meltwater and rain inside the pores has a volume V_MF with unit area, height h_MF, mass M_MF and constant density ${\mathit{\rho }}_i=\mathrm{917}\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$. It is also possible to define the bulk snow density ρ, snow water equivalent SWE and volumetric liquid water content θ_W as $\mathit{\rho }=({\mathit{\rho }}_{\mathrm{S}}{h}_{\mathrm{S}}+{\mathit{\rho }}_{\mathrm{W}}{h}_{\mathrm{W}}+{\mathit{\rho }}_{\mathrm{i}}{h}_{\mathrm{MF}})/h$, $\mathrm{SWE}=\left(\mathit{\rho }h\right)/{\mathit{\rho }}_{\mathrm{W}}$ and ${\mathit{\theta }}_{\mathrm{W}}=h_{\mathrm{W}}/h$, where h is the height of the snowpack equal to $h=h_{\mathrm{S}}+\langle h_{\mathrm{MF}}+h_{\mathrm{W}}-\mathit{\varphi }{h}_{\mathrm{S}}\rangle$ (Avanzi et al., 2015) with ϕ being the porosity and 〈〉 denoting the Macaulay brackets that provide zero when the argument is negative and its value when it is positive. The height h and h_S always coincide except at the end of the snowpack existence when the liquid part and the solid part due to refreezing become predominant (i.e. $h_{\mathrm{MF}}+h_{\mathrm{W}}>\mathit{\varphi }{h}_{\mathrm{S}}$). In this case, h>h_S because a layer of water and/or ice forms on top of the deformable skeleton.

The model solves the mass balance for the dry and liquid mass of the snowpack and the momentum balance and rheological equation for the solid deformable skeleton, resulting in four ordinary differential equations (ODEs) with the variables h_S, h_W, h_MF and ρ_S. The mass fluxes considered are (1) solid precipitation events, snowmelt and wind erosion for the dry-snow mass; (2) rain events, snowmelt, melt–freeze inside the snowpack and run-off for the liquid mass; and (3) melt–freeze for the mass of ice. The dry-snow density is obtained considering (1) the compaction of snow due to compaction not driven by wind, (2) the increase in the densification rate due to drifting snow settlement and (3) densification due to the addition of new mass. The following system is thus obtained (see Appendix Sects. A1–A2 for the derivation of the system and a detailed description of the terms in the equations):

$$\begin{array}{}\text{(1a)}& {\displaystyle \frac{\mathrm{d}{h}_{\mathrm{S}}}{\mathrm{d}t}}=-{\displaystyle \frac{h_{\mathrm{S}}}{{\mathit{\rho }}_{\mathrm{S}}}}{\displaystyle \frac{\mathrm{d}{\mathit{\rho }}_{\mathrm{S}}}{\mathrm{d}t}}+{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{NS}}}{{\mathit{\rho }}_{\mathrm{S}}}}s-(I\cdot a)(T_{\mathrm{A}}-T_{\mathit{\tau }})-{\displaystyle \frac{Q}{{\mathit{\rho }}_{\mathrm{S}}}},\text{(1b)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\mathrm{d}{h}_{\mathrm{W}}}{\mathrm{d}t}}& =r+{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{S}}}{{\mathit{\rho }}_{\mathrm{W}}}}(I\cdot a)(T_{\mathrm{A}}-T_{\mathit{\tau }})+(I^{*}\cdot e\cdot a)(T_{\mathrm{A}}-T_{\mathit{\tau }})\\ & -\mathit{\alpha }\cdot K_{\mathrm{W}},\end{array}\text{(1c)}& {\displaystyle \frac{\mathrm{d}{h}_{\mathrm{MF}}}{\mathrm{d}t}}=-{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{W}}}{{\mathit{\rho }}_i}}(I^{*}\cdot e\cdot a)(T_{\mathrm{A}}-T_{\mathit{\tau }}),\text{(1d)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\mathrm{d}{\mathit{\rho }}_{\mathrm{S}}}{\mathrm{d}t}}& =(c\cdot A_{\mathrm{1}}\cdot U){\mathit{\rho }}_{\mathrm{S}}\exp (-B\cdot (T_{\mathit{\tau }}-T_{\mathrm{S}})-A_{\mathrm{2}}\cdot {\mathit{\rho }}_{\mathrm{S}})\\ & +{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{NS}}-{\mathit{\rho }}_{\mathrm{S}}}{h_{\mathrm{S}}}}s.\end{array}\end{array}$$

In Eq. (1a), ρ_NS is the density of fresh snow (kg m⁻³), s is the solid precipitation rate (m h⁻¹), a is a calibration parameter ($\mathrm{m}{\mathrm{h}}^{-\mathrm{1}}{}^{\circ }{\mathrm{C}}^{-\mathrm{1}}$), T_A and T_τ are the air temperature and the threshold temperature for melting (^∘C), I is equal to $\frac{h_{\mathrm{S}}}{h_{\mathrm{S}}+k}$ with k equal to 0.01 m if T_A≥T_τ and zero otherwise (Avanzi et al., 2015), and Q is the mass of snow eroded by wind ($\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$). In Eq. (1b), r is the liquid precipitation rate (m h⁻¹), e is a calibration parameter, I^* is equal to $\frac{h_{\mathrm{W}}}{h_{\mathrm{W}}+k}$ if T_A<T_τ and to $\frac{h_{\mathrm{MF}}}{h_{\mathrm{MF}}+k}$ if T_A>T_τ (Avanzi et al., 2015), $\mathit{\alpha }=\mathrm{1.9692}\times {\mathrm{10}}^{\mathrm{9}}{\mathrm{m}}^{-\mathrm{1}}{\mathrm{h}}^{-\mathrm{1}}$ (DeWalle and Rango, 2008), and K_W is the intrinsic permeability of water in snow (m²). In Eq. (1d), $c=\mathrm{0.10}\cdot \mathrm{3600}\mathrm{s}{\mathrm{h}}^{-\mathrm{1}}$, A₁=0.0013 m⁻¹, $A_{\mathrm{2}}=\mathrm{0.021}{\mathrm{m}}^{\mathrm{3}}{\mathrm{kg}}^{-\mathrm{1}}$, $B=\mathrm{0.08}{\mathrm{K}}^{-\mathrm{1}}$ (Liston et al., 2007), U is the wind speed contribution (m s⁻¹), and T_S is the average snow temperature (^∘C) obtained assuming thermal equilibrium between the constituents and a bilinear profile of temperature through depth (see De Michele et al., 2013, for further details).

With respect to the model by De Michele et al. (2013) and Avanzi et al. (2015), the version presented in this work includes the contribution of wind erosion to mass balance and effect of wind on densification. This is important when the model is applied to high-altitude sites: Haeberli and Alean (1985), in fact, suggested that a major part of the decrease in accumulation with altitude in the Alps that occurs above about 3500 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ may be due to wind effects.

In analogy with solid transport, snow is mobilized only when wind velocity at the surface exceeds a given threshold that depends on physical properties of the surface snowpack (Li and Pomeroy, 1997). Once transport begins, snow can travel in two main modes: saltation and suspension (Déry and Taylor, 1996; Pomeroy et al., 1997). The total snow transport Q is computed by the model with the following assumptions: (1) only snow erosion occurs, and no deposition of snow eroded in other positions is present; (2) measured wind speed is always referred to a 10 m height – i.e. the height of the snow on the ground is neglected; and (3) wind cannot erode snow that has experienced a temperature greater than 0 ^∘C for the presence of ice crusts or wet layers, following Vionnet et al. (2018). These last two assumptions allow us to compute the series of total snow transport Q decoupled from the snow model since knowledge of snow height is not required. In order to implement the routine, we followed, with some modifications, Lehning et al. (2000), where a model of snowdrift was added to the one-dimensional snow model SNOWPACK (further details about the implementation of the routine are reported in Appendix Sect. A1).

2.2 Model of snow–firn dynamics

We propose here two versions of the snow–firn model. The first version (single-layer) models firn and underlying glacier ice as a single layer (Fig. 1, left panel). The resulting output is an average density and the total column height. The second version (multi-layer) considers separately each firn layer (Fig. 1, right panel), and it allows us to distinguish between layers of firn and glacier ice. The discretized profile of density with depth can be obtained from this second implementation. The snow layer is, instead, treated as a single layer in both versions.

Figure 1 A column of snow, firn and ice as modelled by the single-layer (a) and multi-layer (b) version of the snow–firn model.

In the model we neglected the amount of water percolation inside firn. The presence of water inside firn varies greatly depending on the type of glacier. At high altitudes, where maximum temperatures are rarely positive, the effects of percolation due to melting are limited (Smiraglia et al., 2000); at the cold site of Colle Gnifetti, where the model was applied, percolation occurs only in the few centimetres below the surface and does not involve previous-year layers (Alean et al., 1983). If needed, the structure of the model allows us to easily implement additional processes.

In order to separate snow from firn, we refer to firn's original definition, according to which firn is snow that has survived one melt season (Cuffey and Paterson, 2010).

2.2.1 One-layer modelling of firn

The model is composed of two layers: the snowpack (see Sect. 2.1) and the column of firn below it. The firn is modelled as a single impermeable layer of volume V_F, unit area, height h_F, mass M_F and density ρ_F (Fig. 1, left panel).

The model consists of six ODEs: the four equations of the snow model and in addition the mass balance and momentum balance of firn. The mass variation in firn is obtained considering firn melt, the effects of precipitation on firn and the transformation of snow into firn at the end of each hydrological year. The firn densification rate is obtained considering densification due to overburden stress and densification due to addition of new mass. The resulting system is thus as follows (see Appendix A for the derivation of the system and a detailed description of the terms in the equations):

$$\begin{array}{}\text{(2a)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\mathrm{d}{h}_{\mathrm{S}}}{\mathrm{d}t}}& =-{\displaystyle \frac{h_{\mathrm{S}}}{{\mathit{\rho }}_{\mathrm{S}}}}{\displaystyle \frac{\mathrm{d}{\mathit{\rho }}_{\mathrm{S}}}{\mathrm{d}t}}+{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{NS}}}{{\mathit{\rho }}_{\mathrm{S}}}}s-(I\cdot a)(T_{\mathrm{A}}-T_{\mathit{\tau }})-{\displaystyle \frac{Q}{{\mathit{\rho }}_{\mathrm{S}}}}\\ & -\sum _i{h}_{\mathrm{S}}\mathit{\delta }(t-t_i),\end{array}\text{(2b)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\mathrm{d}{h}_{\mathrm{W}}}{\mathrm{d}t}}& =r+{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{S}}}{{\mathit{\rho }}_{\mathrm{W}}}}(I\cdot a)(T_{\mathrm{A}}-T_{\mathit{\tau }})+(I^{*}\cdot e\cdot a)(T_{\mathrm{A}}-T_{\mathit{\tau }})\\ & -\mathit{\alpha }\cdot K_{\mathrm{W}}-\sum _i{h}_{\mathrm{W}}\mathit{\delta }(t-t_i),\end{array}\text{(2c)}& {\displaystyle \frac{\mathrm{d}{h}_{\mathrm{MF}}}{\mathrm{d}t}}=-{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{W}}}{{\mathit{\rho }}_{\mathrm{i}}}}(I^{*}\cdot e\cdot a)(T_{\mathrm{A}}-T_{\mathit{\tau }})-\sum _i{h}_{\mathrm{MF}}\mathit{\delta }(t-t_i),\text{(2d)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\mathrm{d}{h}_{\mathrm{F}}}{\mathrm{d}t}}& =-{\displaystyle \frac{h_{\mathrm{F}}}{{\mathit{\rho }}_{\mathrm{F}}}}{\displaystyle \frac{\mathrm{d}{\mathit{\rho }}_{\mathrm{F}}}{\mathrm{d}t}}-(I_{\mathrm{F}}\cdot a)(T_{\mathrm{A}}-T_{\mathit{\tau }})\mathit{\delta }\left(h_{\mathrm{S}}\right)\\ & +{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{W}}}{{\mathit{\rho }}_{\mathrm{F}}}}r\mathit{\delta }\left(h_{\mathrm{S}}\right)\langle T_{\mathit{\tau }}-T_{\mathrm{A}}\rangle +\sum _i{\displaystyle \frac{\mathit{\rho }}{{\mathit{\rho }}_{\mathrm{F}}}}h\mathit{\delta }(t-t_i),\end{array}\text{(2e)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\mathrm{d}{\mathit{\rho }}_{\mathrm{S}}}{\mathrm{d}t}}& =(c\cdot A_{\mathrm{1}}\cdot U){\mathit{\rho }}_{\mathrm{S}}\exp (-B\cdot (T_{\mathit{\tau }}-T_{\mathrm{S}})-A_{\mathrm{2}}\cdot {\mathit{\rho }}_{\mathrm{S}})\\ & +{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{NS}}-{\mathit{\rho }}_{\mathrm{S}}}{h_{\mathrm{S}}}}s,\end{array}\text{(2f)}& {\displaystyle \frac{\mathrm{d}{\mathit{\rho }}_{\mathrm{F}}}{\mathrm{d}t}}={\displaystyle \frac{\mathrm{d}{\mathit{\rho }}_{\mathrm{F}}}{\mathrm{d}t}}{|}_{\mathrm{comp}}+\sum _i{\displaystyle \frac{\mathit{\rho }-{\mathit{\rho }}_{\mathrm{F}}}{h_{\mathrm{F}}}}h\mathit{\delta }(t-t_i).\end{array}$$

The last terms in Eqs. (2a)–(2c) move, at the end of each melt season, the remaining snowpack (if present) in the firn layer; t_i is the time instant at the end of hydrological year i and δ(.) is the Dirac delta function. In Eq. (2d), I_F is equal to $\frac{h_{\mathrm{F}}}{h_{\mathrm{F}}+k}$, with k specified above if T_A≥T_τ and equal to zero otherwise. In Eq. (2f), $\frac{\mathrm{d}{\mathit{\rho }}_{\mathrm{F}}}{\mathrm{d}t}{|}_{\mathrm{comp}}$ is the densification of firn due to compaction (see Sect. 2.4). Equations (2a)–(2f) are impulsive differential equations (see e.g. Bainov and Simeonov, 1993, for mathematical details). This type of differential equation involving the impulse effect is used to describe the evolution of many physical phenomena that have a sudden change in their states such as mechanical systems with impact; biological systems such as heartbeats and blood flows; and population dynamics.

2.3 Multi-layer modelling of firn

Firn is modelled as a multi-layer column where each layer j has volume V_Fj, unit area, height h_Fj, mass M_Fj and density ρ_Fj.

The equations of the model change as follows:

$$\begin{array}{}\text{(3a)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\mathrm{d}{h}_{\mathrm{S}}}{\mathrm{d}t}}& =-{\displaystyle \frac{h_{\mathrm{S}}}{{\mathit{\rho }}_{\mathrm{S}}}}{\displaystyle \frac{\mathrm{d}{\mathit{\rho }}_{\mathrm{S}}}{\mathrm{d}t}}+{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{NS}}}{{\mathit{\rho }}_{\mathrm{S}}}}s-(I\cdot a)(T_{\mathrm{A}}-T_{\mathit{\tau }})\\ & -{\displaystyle \frac{Q}{{\mathit{\rho }}_{\mathrm{S}}}}-\sum _i{h}_{\mathrm{S}}\mathit{\delta }(t-t_i),\end{array}\text{(3b)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\mathrm{d}{h}_{\mathrm{W}}}{\mathrm{d}t}}& =r+{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{S}}}{{\mathit{\rho }}_{\mathrm{W}}}}(I\cdot a)(T_{\mathrm{A}}-T_{\mathit{\tau }})+(I^{*}\cdot e\cdot a)(T_{\mathrm{A}}-T_{\mathit{\tau }})\\ & -\mathit{\alpha }\cdot K_{\mathrm{W}}-\sum _i{h}_{\mathrm{W}}\mathit{\delta }(t-t_i),\end{array}\text{(3c)}& {\displaystyle \frac{\mathrm{d}{h}_{\mathrm{MF}}}{\mathrm{d}t}}=-{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{W}}}{{\mathit{\rho }}_{\mathrm{i}}}}(I^{*}\cdot e\cdot a)(T_{\mathrm{A}}-T_{\mathit{\tau }})-\sum _i{h}_{\mathrm{MF}}\mathit{\delta }(t-t_i),\text{(3d)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\mathrm{d}{h}_{\mathrm{F}\mathrm{1}}}{\mathrm{d}t}}& =-{\displaystyle \frac{h_{\mathrm{F}\mathrm{1}}}{{\mathit{\rho }}_{\mathrm{F}\mathrm{1}}}}{\displaystyle \frac{\mathrm{d}{\mathit{\rho }}_{\mathrm{F}\mathrm{1}}}{\mathrm{d}t}}-(I_{\mathrm{F}}\cdot a)(T_{\mathrm{A}}-T_{\mathit{\tau }})\mathit{\delta }\left(h_{\mathrm{S}}\right)\\ & +{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{W}}}{{\mathit{\rho }}_{\mathrm{F}\mathrm{1}}}}r\mathit{\delta }\left(h_{\mathrm{S}}\right)\langle T_{\mathit{\tau }}-T_{\mathrm{A}}\rangle +\sum _i{\displaystyle \frac{\mathit{\rho }}{{\mathit{\rho }}_{\mathrm{F}\mathrm{1}}}}h\mathit{\delta }(t-t_i)\\ & -\sum _i{h}_{\mathrm{F}\mathrm{1}}\mathit{\delta }(t-t_i),\end{array}\text{(3e)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\mathrm{d}{h}_{\mathrm{F}j}}{\mathrm{d}t}}& =-{\displaystyle \frac{h_{\mathrm{F}j}}{{\mathit{\rho }}_{\mathrm{F}j}}}{\displaystyle \frac{\mathrm{d}{\mathit{\rho }}_{\mathrm{F}j}}{\mathrm{d}t}}+\sum _i{h}_{{\mathrm{F}}_{j-\mathrm{1}}}\mathit{\delta }(t-t_i)\\ & -\sum _i{h}_{\mathrm{F}j}\mathit{\delta }(t-t_i),\end{array}\text{(3f)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\mathrm{d}{\mathit{\rho }}_{\mathrm{S}}}{\mathrm{d}t}}& =(c\cdot A_{\mathrm{1}}\cdot U){\mathit{\rho }}_{\mathrm{S}}\exp (-B\cdot (T_{\mathit{\tau }}-T_{\mathrm{S}})-A_{\mathrm{2}}\cdot {\mathit{\rho }}_{\mathrm{S}})\\ & +{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{NS}}-{\mathit{\rho }}_{\mathrm{S}}}{h_{\mathrm{S}}}}s,\end{array}\text{(3g)}& {\displaystyle \frac{\mathrm{d}{\mathit{\rho }}_{\mathrm{F}\mathrm{1}}}{\mathrm{d}t}}={\displaystyle \frac{\mathrm{d}{\mathit{\rho }}_{\mathrm{F}\mathrm{1}}}{\mathrm{d}t}}{|}_{\mathrm{comp}}+\sum _i{\displaystyle \frac{\mathit{\rho }-{\mathit{\rho }}_{\mathrm{F}\mathrm{1}}}{h_{\mathrm{F}\mathrm{1}}}}h\mathit{\delta }(t-t_i),\text{(3h)}& {\displaystyle \frac{\mathrm{d}{\mathit{\rho }}_{\mathrm{F}j}}{\mathrm{d}t}}={\displaystyle \frac{\mathrm{d}{\mathit{\rho }}_{\mathrm{F}j}}{\mathrm{d}t}}{|}_{\mathrm{comp}},\end{array}$$

where j goes from 2 to the total number of firn layers. Firn layers that reach the ice density or whose height goes to zero are removed from the model.

2.4 Firn densification

The densification of firn due to compaction is usually subdivided into three stages: (1) a first stage dominated by the settling of grains that allows us to reach densities of up to about 550 kg m⁻³, (2) a second stage dominated by sintering that extends up to the close-off density (i.e. the density at which pores become isolated) of about 830 kg m⁻³ and (3) a last stage that ends when ice density is reached in which further densification is driven by the compression of the bubbles of air (Cuffey and Paterson, 2010). This last stage is in turn subdivided into two phases, depending on if the pores are cylindrical or spherical. Different models of firn densification are available in the literature (see Lundin et al., 2017, for a review). Here we implemented the model of Arnaud et al. (2000) with some of the modifications proposed by Bréant et al. (2017) (we will refer to it with the abbreviation AR) and the model of Herron and Langway (1980) (we will refer to it with the abbreviation HL). Other models could also be implemented. Both HL and AR were developed for polar sites. The HL model was derived using ice cores with a mean annual firn temperature between −57 and −15 ^∘C and a mean annual accumulation between 0.022×10³ and 0.5×10³ $\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$, while the AR model was derived from cores with a mean annual firn temperature between −57 and −19 ^∘C and a mean annual accumulation between 0.022×10³ and 1.1×10³ $\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$. In the model of AR, densification equations are based on grain sliding and creep deformations, even though they maintain empirical parameters. The model of HL consists of empirical equations tuned with ice cores, based on the assumption that a proportionality is present between the variation in density and the variation in stress due to new accumulation. Besides, the model of AR represents stresses explicitly, while in HL the load is parametrized through annual surface accumulation.

The model of HL was already applied to non-polar ice cores by Huss (2013), where the model was recalibrated in order to match depth–density profiles of temperate and polythermal firn. In the presented application, the parameters were not recalibrated, despite the fact that the study site is an alpine site. This was motivated by the low mean annual firn temperature (MAFT) and low surface accumulation observed at Colle Gnifetti that resemble conditions of some polar sites.

In AR all three stages of firn densification are modelled. Equations are as follows:

$$\begin{array}{}\text{(4)}& \begin{array}{rl}& {\displaystyle \frac{\mathrm{d}{\mathit{\rho }}_{\mathrm{F}}}{\mathrm{d}t}}{|}_{\mathrm{comp}}=\\ & \left\{\begin{array}{ll}\mathit{\gamma }\frac{max(P,{\mathrm{10}}^{\mathrm{4}})}{({\mathit{\rho }}_{\mathrm{F}}/{\mathit{\rho }}_i{)}^{\mathrm{2}}}\left(\mathrm{1}+\frac{\mathrm{0.5}}{\mathrm{6}}-\frac{\mathrm{5}}{\mathrm{3}}\frac{{\mathit{\rho }}_{\mathrm{F}}}{{\mathit{\rho }}_i}\right){\mathit{\rho }}_i& D_{\mathrm{D}}\le {\mathit{\rho }}_{\mathrm{F}}/{\mathit{\rho }}_i\le D_{\mathrm{0}}\\ \mathrm{5.3}A\cdot {\left(({\mathit{\rho }}_{\mathrm{F}}/{\mathit{\rho }}_i{)}^{\mathrm{2}}{D}_{\mathrm{0}}\right)}^{\mathrm{1}/\mathrm{3}}{\left(\frac{a_{\mathrm{c}}}{\mathit{\pi }}\right)}^{\mathrm{1}/\mathrm{2}}{\left(\frac{\mathrm{4}\mathit{\pi }\cdot P\cdot {\mathit{\rho }}_i}{\mathrm{3}{a}_{\mathrm{c}}\cdot Z\cdot {\mathit{\rho }}_{\mathrm{F}}}\right)}^{\mathrm{3}}{\mathit{\rho }}_i& D_{\mathrm{0}}<{\mathit{\rho }}_{\mathrm{F}}/{\mathit{\rho }}_i\le D_{\mathrm{c}}\\ \mathrm{2}A\cdot \frac{{\mathit{\rho }}_{\mathrm{F}}(\mathrm{1}-{\mathit{\rho }}_{\mathrm{F}}/{\mathit{\rho }}_i)}{{\mathit{\rho }}_i(\mathrm{1}-(\mathrm{1}-{\mathit{\rho }}_{\mathrm{F}}/{\mathit{\rho }}_i{)}^{\frac{\mathrm{1}}{\mathrm{3}}}{)}^{\mathrm{3}}}{\left(\frac{\mathrm{2}(P-P_{\mathrm{b}})}{\mathrm{3}}\right)}^{\mathrm{3}}{\mathit{\rho }}_i& D_{\mathrm{c}}<{\mathit{\rho }}_{\mathrm{F}}/{\mathit{\rho }}_i\le \mathrm{0.95}\\ \frac{\mathrm{9}}{\mathrm{4}}A\cdot (\mathrm{1}-{\mathit{\rho }}_{\mathrm{F}}/{\mathit{\rho }}_i)(P-P_{\mathrm{b}}){\mathit{\rho }}_i& {\mathit{\rho }}_{\mathrm{F}}/{\mathit{\rho }}_i>\mathrm{0.95}.\end{array}\right.\end{array}\end{array}$$

In the first stage ($D_{\mathrm{D}}\le {\mathit{\rho }}_{\mathrm{F}}/{\mathit{\rho }}_i\le D_{\mathrm{0}}$), P is the overburden pressure (Pa) and $\mathit{\gamma }={\mathit{\gamma }}^{\prime }\exp \left(-\frac{Q_{\mathrm{1}}}{R_{\mathrm{G}}(T_{\mathrm{F}}+\mathrm{273.15})}\right)$ in which R_G is the gas constant, Q₁ an activation energy equal to $\mathrm{48}\times {\mathrm{10}}^{\mathrm{3}}\mathrm{J}{\mathrm{mol}}^{-\mathrm{1}}$, T_F the average temperature of firn (^∘C), and γ^′ a parameter whose value is set in order to have a continuous densification rate between the first and second stage (estimated in Sect. 4.2). D_D is the relative surface density, and D₀ is the relative density at the transition between the first stage and the second stage. In the second stage ($D_{\mathrm{0}}<{\mathit{\rho }}_{\mathrm{F}}/{\mathit{\rho }}_i\le D_{\mathrm{c}}$), $A=A_{\mathrm{0}}\exp \left(-\frac{Q_{\mathrm{2}}}{R_{\mathrm{G}}(T_{\mathrm{F}}+\mathrm{273.15})}\right)$ with $A_{\mathrm{0}}=\mathrm{2.84}\times {\mathrm{10}}^{-\mathrm{11}}$ ${\mathrm{Pa}}^{-\mathrm{3}}{\mathrm{h}}^{-\mathrm{1}}$, a_c is the average contact area, Z is the number of particle contacts (see Appendix Sect. A3 for the expression of a_c and Z) and Q₂ is an activation energy. The value of Q₂ was set to 60×10³ J mol⁻¹, as in the model of Arnaud et al. (2000), since it is the typical activation energy associated with self-diffusion of ice. However, at higher temperature (i.e. higher than −10 ^∘C) a higher activation energy may be required to best fit density profiles with firn densification models (Cuffey and Paterson, 2010; Arthern et al., 2010; Jacka and Jun, 1994). A discussion of the thermal variation in the creep parameter and the impact of the different sintering mechanisms on it can be found in Bréant et al. (2017). Lastly, in the third stage (${\mathit{\rho }}_{\mathrm{F}}/{\mathit{\rho }}_i>D_{\mathrm{c}}$), P_b is the pressure inside the bubbles equal to $P_{\mathrm{b}}=P_{\mathrm{c}}\frac{({\mathit{\rho }}_{\mathrm{F}}/{\mathit{\rho }}_i)(\mathrm{1}-D_{\mathrm{c}})}{D_{\mathrm{c}}\cdot (\mathrm{1}-{\mathit{\rho }}_{\mathrm{F}}/{\mathit{\rho }}_i)}$ with D_c and P_c the relative density and pressure at the transition between the second and third stage. In the single-layer version, the overburden pressure P was computed as the overburden of the snowpack layer plus half of the firn layer. In the multi-layer version, we computed the overburden for each layer of firn as the overburden of the snowpack plus the overburden of all the firn layers above plus the overburden of half the firn layer considered.

In HL only the first and second densification stages are modelled. The equations are as follows:

$$\begin{array}{}\text{(5)}& \begin{array}{rl}& {\displaystyle \frac{\mathrm{d}{\mathit{\rho }}_{\mathrm{F}}}{\mathrm{d}t}}{|}_{\mathrm{comp}}=\\ & \left\{\begin{array}{ll}{k}_{\mathrm{0}}\cdot (\mathit{\omega }\times {\mathrm{10}}^{\mathrm{3}})\cdot ({\mathit{\rho }}_i-{\mathit{\rho }}_{\mathrm{F}})& {\mathit{\rho }}_{\mathrm{D}}\le {\mathit{\rho }}_{\mathrm{F}}\le \mathrm{550}\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}\\ k_{\mathrm{1}}\cdot (\mathit{\omega }\times {\mathrm{10}}^{\mathrm{3}}{)}^{\mathrm{0.5}}\cdot ({\mathit{\rho }}_i-{\mathit{\rho }}_{\mathrm{F}})& \mathrm{550}<{\mathit{\rho }}_{\mathrm{F}}<\mathrm{800}\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}\end{array}\right.\end{array}\end{array}$$

where $k_{\mathrm{0}}=\mathrm{11}\exp \left(-\frac{\mathrm{10}\mathrm{160}}{R_{\mathrm{G}}(T_{\mathrm{F}}+\mathrm{273.15})}\right)$, $k_{\mathrm{1}}=\mathrm{575}\exp \left(-\frac{\mathrm{21}\mathrm{400}}{R_{\mathrm{G}}(T_{\mathrm{F}}+\mathrm{273.15})}\right)$ and ω is the annual snow accumulation ($\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$). In HL the transition density between the first and second stage is fixed and equal to 550 kg m⁻³. In order to run the model of HL in a dynamic way, for each year we computed the annual accumulation averaging the ones modelled between the year of deposition of the firn layer and the year before the one considered, following Stevens et al. (2020).

Steady-state firn densification models are not applied to the superficial snow where the metamorphism is more complex and significantly influenced by air temperature. The original model of Arnaud et al. (2000), for example, was used only for depths higher than 2 m. In this case, we applied them only to densities higher than a density ρ_D, which represents the average snow density. For firn densities lower than ρ_D, the densification equation of snow was adopted although neglecting wind contribution. In this way, the transition between the two equations is driven by density rather than associated with the end of a water year. This is important, for example, when consistent fresh snow falls over the snowpack at the end of the hydrological year.

2.4.1 Temperature profile

The energetic description of the volume was simplified assuming the constituents were in thermal equilibrium and assuming a bilinear profile of temperature through depth. Temperature was assumed to vary linearly from the surface temperature T₀ to the MAFT at the depth z_M at which seasonal variation in temperature is negligible. Below z_M, temperature was kept constant and equal to MAFT. In cold glaciers the value of MAFT is close to the mean annual air temperature (MAAT) when meltwater percolation is limited (Suter et al., 2001) while in temperate glaciers it is equal to the melting temperature (Cuffey and Paterson, 2010). Surface temperature was fixed equal to T_A if T_A<0 ^∘C and zero otherwise. Huss (2013) has already assumed a bilinear profile of temperature in order to study temperate firn densification, fixing z_M to 5 m since it is the typical penetration of winter air temperature. The temperature profile was then used to compute the average snow and firn temperatures that influence snow and firn densification.

2.5 Numerical model

The model was solved using the forward Euler method with a constant step size, Δt, of 1 h. To also compute the last terms in Eqs. (1d), (2f) and (3g) when h_S, h_F and $h_{{\mathrm{F}}_{\mathrm{1}}}$ are zero, these terms were calculated, following De Michele et al. (2013), as $\frac{{\mathit{\rho }}_{\mathrm{NS}}\left(t\right)-{\mathit{\rho }}_{\mathrm{S}}}{h_{\mathrm{S}}\left(t\right)+s\left(t\right)\mathrm{\Delta }t}s\left(t\right)$, $\frac{\mathit{\rho }\left(t\right)-{\mathit{\rho }}_{\mathrm{F}}}{h_{\mathrm{F}}\left(t\right)+h\left(t\right)}\frac{h\left(t\right)}{\mathrm{\Delta }t}$ and $\frac{\mathit{\rho }\left(t\right)-{\mathit{\rho }}_{{\mathrm{F}}_{\mathrm{1}}}}{h_{{\mathrm{F}}_{\mathrm{1}}}\left(t\right)+h\left(t\right)}\frac{h\left(t\right)}{\mathrm{\Delta }t}$. Regarding the vertical discretization, the firn component of the multi-layer version of the model was discretized modelling one layer for each hydrological year.

3 Study area and data

In the following section we will present the study area (Sect. 3.1), the data collection and handling (Sect. 3.2–3.3), and finally the calibration and site-specific parameters (Sect. 3.4–3.5).

3.1 Study area

The site of Colle Gnifetti (CG) is part of the summit ranges of the Monte Rosa massif, Swiss–Italian Alps. It is the uppermost part of the accumulation area of the Grenzgletscher (Border Glacier), and it forms a saddle that lies between the Signalkuppe (4554 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$) and Zumsteinspitze (4563 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$) at an altitude of 4400–4550 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ (Lüthi and Funk, 2000) (Fig. 2). The glacier at Colle Gnifetti has a thickness of between 60 and 120 m and a MAFT of −14 ^∘C (Wagenbach et al., 2012). The regime is that of a high-altitude site, i.e. nearly persistent sub-zero air temperature, a high precipitation total and high wind speed (Suter et al., 2001). A mean annual precipitation of 2.7×10³ $\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ with an interannual variability of $\mathrm{0.8}\times {\mathrm{10}}^{\mathrm{3}}\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ (Mariani et al., 2014) was estimated for the period 1961–1993 from a core extracted at the upper Grenzgletscher (Eichler et al., 2000).

Figure 2 The site of Colle Gnifetti and the location of the ice cores considered in the present work. CG03 and CG15 ice core share the same location therefore CG03 is not shown. The position of Capanna Regina Margherita (CM) is also shown.

Even though the site is characterized by high precipitation totals, accumulation in the saddle is considerably lower and highly variable over the glacier surface due to wind erosion, with values ranging from about 0.15×10³ to 1.2×10³ $\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ depending on the wind exposure (Alean et al., 1983; Lüthi and Funk, 2000; Licciulli et al., 2020). Alean et al. (1983) measured the accumulation at CG between 17 August 1980 and 23 July 1982 with a network of 30 stakes. For the period between 14 August 1981 and 23 July 1982 the mass balance was negative at all the stakes due to wind erosion, while the net accumulation of the hydrological year 1980–1981 varied between $+\mathrm{0.04}\times {\mathrm{10}}^{\mathrm{3}}$ and $+\mathrm{1.18}\times {\mathrm{10}}^{\mathrm{3}}$ $\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ with the highest values on south-facing slopes. This occurs because the enhanced melting and refreezing cause the formation of wet layers and ice crusts and because higher temperatures are associated with faster densification, and both these aspects reduce the possibility of wind eroding snow. This also results in the fact that almost all the snow that survives the melt season comes from summer events (Bohleber et al., 2018; Schöner et al., 2002).

3.2 Data collection

The stations whose data were used in this study are presented in Fig. 3, and they are summarized in Table 1. Hourly data of air temperature and wind speed at Capanna Regina Margherita (CM, Margherita Hut in English) were used as input for the model; hourly data at Passo del Moro (PM, Monte Moro Pass in English) to reconstruct precipitation at CM; and hourly and daily air temperature data at Macugnaga Pecetto (MP), Passo del Moro, Bocchetta delle Pisse (BDP) and Ceppo Morelli (CPM) to infill missing temperature data at Capanna Regina Margherita. Hourly wind speed data at Valtournenche–Cime Bianche (CB) and Col du Grand St Bernard (SB, Great St Bernard Pass in English) were used to infill missing wind speed data at CM. Hourly data at the meteorological stations of Gressoney-Saint-Jean – Weissmatten (GWm) and Gressoney-la-Trinité – Gabiet (GGm) and snow water equivalent data (GWs and GGs) were used to calibrate and validate the parameters a and e of the snow model. Snow water equivalent is measured by the Aosta Valley Region during winter both at fixed locations and at itinerant sites. For GGs, 4 years of measurements was available with on average 24 data points for each winter. For GWs, 5 years was available with on average 8 data points for each winter.

The station of Capanna Regina Margherita, whose data were used to run the snow–firn model, was installed in 2002 by the Piedmont Region at the Regina Margherita Hut as part of a project that aimed to study the interaction between synoptic flow and orography. With its 4560 m of altitude, it can be considered the highest meteorological station in Europe, and its wind speed series can be considered representative of the synoptic conditions (Martorina et al., 2003). Due to its recent installation, the use of these data limits the length of the simulation and the number of cores with which our results can be compared. Nevertheless, we believe that, given the peculiar characteristics of the station, the use of these data may give added value to this study.

Figure 3 Location of the meteorological stations used: Capanna Regina Margherita (CM), Macugnaga Pecetto (MP), Ceppo Morelli (CPM), Passo del Moro (PM), Bocchetta delle Pisse (BDP), Col du Grand St Bernard (SB), Valtournenche–Cime Bianche (CB), Gressoney-la-Trinité – Gabiet (GGm and GGs) and Gressoney-Saint-Jean – Weissmatten (GWm and GWs). The location of the meteorological station at Weissmatten and the location of snow measurements are only a few metres apart, so only one point is reported (identified with GW). In the bottom left panel a zoom over some stations is included. All stations belong to Arpa Piemonte with the exclusion of CB, GGm, GGs, GWm and GWs, which belong to the Aosta Valley Region, and SB, which belongs to the National Oceanic and Atmospheric Administration (NOAA). Source of the basemap: Arpa Piemonte geoportal.

In Table 2 ice core data are reported (Fig. 2). Available data consist of some or all of the following information: depth in metres, depth in metres of water equivalent, density and dating. We recall that the first three variables are related, so one of them can be computed given the other two.

Table 1 Meteorological data employed in the case study (p stands for precipitation, SD for snow depth, T_A for air temperature, u for average wind speed and s for fresh snow). Hydrological years are identified by the last year; e.g. 2009 is hydrological year 2008–2009. With the term hydrological year we refer to the period from 1 October to 30 September of the next year.

Gäggeler et al. (1983)Gäggeler et al. (1983)Sigl et al. (2018)Sigl et al. (2018)Ardenghi (2012)Licciulli et al. (2020)Licciulli et al. (2020)Licciulli et al. (2020)

Table 2 Ice core data employed in the case study.

3.3 Data handling

The model requires as input a continuous series of air temperature, precipitation and wind speed.

Following the comparison presented by Henn et al. (2013), to fill missing hourly temperature data at Capanna Margherita, the MicroMet preprocessor (Liston and Elder, 2006) was adopted for gaps smaller than 24 h and a long-term lapse rate approach with five stations (CM, MP, CPM, PM, BDP) was adopted for longer gaps. MicroMet is a meteorological model that includes a data-fill procedure adopted here. The method distinguishes between three conditions: (1) for 1 h gaps the missing information is replaced with the average of the previous and next measurement; (2) for 2–24 h gaps each missing value is replaced with the average of the values recorded the next and previous day at the same hour; (3) for longer gaps an auto-regressive integrated moving average (ARIMA) model is used (Hyndman and Athanasopoulos, 2021). Here, for the third condition, the long-term lapse rate approach was used, as specified above. In the period 1 October 2002–13 August 2013, 0.37 % of hourly temperature data were missing. After the MicroMet procedure 0.23 % remained missing and were substituted with a long-term lapse rate approach.

Wind speed data measured at CM are characterized by repeated zero values that are not observed in nearby stations and that are probably due to the freezing of the anemometer. In the period 1 October 2002–13 August 2013 nearly 30 % of the wind speed data at CM were equal to zero, while 1.3 % were missing. By comparison, in the same period, there were 2 % zero values in SB series. These zero values were therefore considered missing. To fill missing wind speed data at CM, the MicroMet procedure was used for gaps smaller than 24 h. For gaps longer than 24 h, data were replaced using measurements at CB or, if wind speed data were also missing at CB, with data measured at SB. In both series, zero wind speed values recorded for more than 4 consecutive hours were set as missing. In order to take into account the different characteristics of the sites, we first computed for each of the three stations the mean and standard deviation for each hour of the year, and we removed this from the data. Missing data at CM were first replaced with the corresponding residual value measured at CB (or SB), and then the final value was obtained using the mean and standard deviation estimated at CM. Reconstructed negative wind speed values at CM were set to zero. Missing wind speed data at CM were set to zero if they were zero at CB (or SB).

The precipitation series at CM was reconstructed using hourly data measured at PM. The station was chosen due to it being in the vicinity of CM and its altitude of 2820 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}$. Using the formula proposed by Alpert (1986) and considering a bell-shaped mountain, we estimated for the Monte Rosa massif an altitude of maximum precipitation of z_m=2547 m. The altitude of maximum precipitation is away from the crest as is typical of large mountains (Roe, 2005). We therefore expect to have similar precipitation totals at CM and PM. The precipitation series measured at PM needs to be integrated with snow depth data due to the under-catch of solid precipitation by the pluviometer or does not catch solid precipitation events in winter. In order to reconstruct the total precipitation series, we followed the routine presented by Avanzi et al. (2014). Solid precipitation is obtained looking at the positive variations in snow depth data, while rainfall is given by the difference, if positive, between total precipitation and solid precipitation. Positive variations in snow depth, however, may also be recorded when strong temperature variations occur, thus introducing false events. Unlike Avanzi et al. (2014), we approached this problem smoothing the snow depth series with a moving average whose window size was calibrated running the snow model at PM and looking for the best match between simulations and observations. Even though PM has an altitude higher than the estimated altitude of maximum precipitation, we obtained a mean annual precipitation of about 2×10³ $\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ for the period 2002–2019, lower than the one estimated at CG by Mariani et al. (2014). We suppose this may be due to wind erosion events; the procedure implemented by Avanzi et al. (2014), in fact, may compensate for snow depth variations due to wind erosion decreasing the estimated solid precipitation. We therefore increased the resulting hourly solid precipitation with a constant factor in order to match the observed mean annual accumulation at CG of 2.7×10³ $\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$. The total precipitation series was then divided between solid and liquid precipitation using a threshold of 1 ^∘C since this is the value generally found in Europe (Jennings et al., 2018).

3.4 Calibration of model's parameters

The model requires the calibration of three parameters, namely a and e in Eqs. (1a)–(1c) and γ^′ in Eq. (4).

The parameters a and e were calibrated running the snow model, with the addition of the wind module, at Gressoney-Saint-Jean – Weissmatten, with an hourly time step from 1 October 2015 to 30 September 2020. Input series were processed as reported for PM in Sect. 3.3. The parameter a governs the amount of snowmelt and consequently snow height and the relative amount of snow and ice inside the snowpack, thus influencing snow water equivalent and density. On the contrary, the parameter e influences only the relative amount of snow and ice inside the snowpack and does not contribute to snow height. The calibration problem is therefore a multi-objective one, since we could optimize the error in both snow height and density or SWE. We decided to move from a multi-objective to a single-objective optimization problem, aggregating the NSE (Nash–Sutcliffe efficiency) between observed and simulated snow depth data and the NSE between observed and simulated snow water equivalent data. We calculated the error metrics considering together all the available years but computing the measure only in the periods with snow depth higher than zero, and we aggregated them, giving a weight of 0.7 to the first and 0.3 to the second. In this way we took into account the higher uncertainty in SWE data due to the shorter sample length and the non-coincidence between the location of the meteorological station and the snow measurements. The optimum parameters were then estimated for different moving-average windows, used to process solid precipitation input data, with the use of a population-evolution-based algorithm, namely SCE-UA (Shuffled Complex Evolution – University of Arizona) (Duan et al., 1992, 1993). We thus obtained a pair of a and e values for each window, and we selected the one maximizing the objective function. The validation was performed applying the model with the selected set of parameters at Gressoney-la-Trinité – Gabiet for the period 1 October 2017–30 September 2021.

The parameter γ^′, which governs the firn densification rate in AR, was chosen in order to have a continuous densification rate between the first and second stage of densification. For each of the available ice cores, with the exception of CG11, we computed the parameter γ^′ running AR in a steady-state condition (Bader, 1954) using the mean accumulation reported in Table 2. In addition, the parameters of the firn densification model chosen may need calibration if applied to sites significantly different from the polar ones.

3.5 Site-specific parameters

In order to apply AR we use D₀=0.56 (Bréant et al., 2017), $P_{\mathrm{c}}=\mathrm{740}\times {\mathrm{10}}^{\mathrm{2}}\mathrm{Pa}$ (Lüthi and Funk, 2000) and D_c=0.9 since the precise value is not known at CG (Lüthi and Funk, 2000). Two different firn temperatures, $T_{\mathrm{F}}=-\mathrm{14}{}^{\circ }\mathrm{C}$ and $T_{\mathrm{F}}=-\mathrm{10}{}^{\circ }\mathrm{C}$, that cover the observed ice temperatures at CG (Lüthi and Funk, 2000) were tested together with different surface densities, chosen looking at values already used in the literature at CG. We selected three values: ρ_D=300 kg m⁻³, ρ_D=360 kg m⁻³ and ρ_D=410 kg m⁻³, the values already assumed by Licciulli et al. (2020) and Lüthi and Funk (2000). The model of AR was run with a slight modification. We used the first-stage densification equation up to a relative density of 0.6, but we kept D₀=0.56 in the second-stage densification equation. The latter, in fact, cannot be applied for D=D₀, and it gives densification rates tending to infinity for values tending to D₀. The other site-specific parameters of the snow–firn model that require specification are z_M, set to 5 m (Haeberli and Funk, 1991), and the grain radius R, which influences the threshold wind speed. It is defined as $R=\mathrm{3}/\left({\mathit{\rho }}_i\mathrm{SSA}\right)$, where SSA is the specific surface area in m² kg⁻¹. SSA was computed adopting the parametrization of Domine et al. (2007) for recent snow: $\mathrm{SSA}=-\mathrm{16.051}\ln ({\mathit{\rho }}_{\mathrm{S}}\times {\mathrm{10}}^{-\mathrm{3}})+\mathrm{7.01}$.

4 Results

4.1 Parameters' estimation

We obtained a value of the parameters a and e of 2.94×10⁻⁴ m h^−1 ∘C⁻¹ and 0.164, respectively. The combined NSE in calibration is 0.82, with a NSE of 0.84 and 0.78 for snow depth and snow water equivalent data, respectively. In validation, the combined NSE and the NSE values for snow depth and snow water equivalent are 0.77, 0.84 and 0.61, respectively. We also computed the RMSE and mean bias error (MBE) in validation, which are equal to 0.126×10³ and 0.0116×10³ $\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ for snow water equivalent and 0.26 and 0.0672 m for snow depth. Avanzi et al. (2014) estimated the parameter a for a selection of 40 sites with altitudes between 91 and 3389 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ within the SNOTEL (Snow Telemetry) network, a network of automated stations located in mountain basins of the western USA and operated by the Natural Resources Conservation Service (NRCS). They obtained median values of a of between 1×10⁻⁴ and 6×10⁻⁴ m h^−1 ∘C⁻¹. Regarding the parameter e, values of 0.2 and 0.25 were estimated by Avanzi et al. (2015) for two sites in Japan.

4.2 Steady-state firn densification

The depth–density profiles obtained using the model of AR and HL in a steady-state condition are reported in Fig. 4 for a surface density ρ_D=360 kg m⁻³. Both HL and AR have very good performances when applied to the KCC ice core. The worst performances occur for CG03 with an underestimation of the densification rate for all depths. For the remaining ice cores the models of AR and HL have a good fit up to depths of about 20–30 m, but they underestimate the densification rate below it. The profiles show in general a better performance of HL. We recall that the model of HL was derived also considering cores with MAFT and accumulation close to the ones of CG, while AR was optimized for cores with lower MAFT.

Figure 4 Observed and modelled depth–density profiles. Modelled profiles are obtained running Arnaud et al. (2000) (AR) and Herron and Langway (1980) (HL) in a steady-state condition. The numbers 14 and 10 stand for a mean annual firn temperature of −14 and −10 ^∘C, respectively.

4.3 Snow accumulation

The annual accumulation obtained from the snow–firn model is reported in Fig. 5, along with the values retrieved from the three available ice cores, the average value of the observations and its 95 % confidence interval. The RMSE between the model and the average of the observations is equal to 0.22×10³ $\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$, and the modelled and observed average annual accumulation and standard deviation are reported in Table 3.

Figure 5 Annual accumulation modelled and retrieved from three ice cores. The average of the annual accumulations from ice cores and its 95 % confidence interval are also reported.

Table 3 Modelled and observed mean (μ) and standard deviation (σ) of the accumulation rate for the period 2003–2010.

In order to better understand the characteristics of the accumulation at CG, the monthly box plots of snow transport; solid precipitation; number of hours with $T_{\mathrm{A}}>\mathrm{0}{}^{\circ }\mathrm{C}$, which in the model correspond to hours with melting; and monthly contribution to annual accumulation, computed for each month as 100 × (solid precipitation − eroded snow) $/$ solid precipitation, are provided in Fig. 6. Since snow is moved into firn at the end of September and wind is not allowed to erode firn, the fraction of conserved snow of September may be overestimated and the snow transport of October underestimated. We can see that annual accumulation is composed of snow deposited mainly between April and September, with June the month that on average contributes the most. The months in which solid precipitation is conserved are also the months in which temperature goes above the melting point; winter snow, instead, is completely removed.

Figure 6 Box plots of monthly eroded snow (a), monthly solid precipitation (b), monthly number of hours with above-zero temperatures (c) and monthly fraction of conserved solid precipitation (d), obtained for the period 1 October 2002–30 September 2015. The horizontal bar inside the box is drawn at the median, and the upper and lower ends of the box are drawn at the upper and lower quartile, respectively. The vertical lines, called whiskers, extend up to the most distant point that has a value within 1.5 times the interquartile range. The points outside these limits are drawn individually with dots.

4.4 Firn density

4.4.1 One-layer model version

The modelled firn density was compared with the density estimated from KCC and CG15 ice cores. With the one-layer version, we obtain one average value of firn density for each run of the model. We therefore run the model multiple times, fixing the end year of the simulation to the date of the core drilling and anticipating at each run the starting date of 1 year. For each run, the corresponding observed firn density was obtained averaging the density profile of the ice core associated with the same range of years. The results obtained implementing both AR and HL are reported in Fig. 7. For KCC we fixed the MAFT to −14 ^∘C, while for CG15 we fixed it to −10 ^∘C, looking for the best performance in Fig. 4.

Regarding firn density, we have contrasting results depending on the core and the densification model adopted. Both model versions overestimate KCC density with a better performance when AR is implemented; on the contrary we obtained an underestimation of CG15 average density, with a better performance when HL is implemented. In all the combinations, however, we observed a reduction in the error when more firn layers are averaged.

Moving to firn depth, the model nearly always predicts higher depths, with more significant differences for the KCC ice core.

Figure 7 Observed and modelled (with one-layer version) average firn density (dots, right axis) and depth (bars, left axis) for KCC and CG15. Each cluster corresponds to a run of the model whose starting date is reported on the x axis. The ending year for all runs is 30 August 2013 for KCC and 30 September 2015 for CG15. AR and HL stand for the model versions with Arnaud et al. (2000) and Herron and Langway (1980) implemented, respectively. The values 300, 360 and 410 stand for the chosen value of ρ_D.

4.4.2 Multi-layer model version

The modelled density profile was compared with KCC and CG15 density data, implementing in the model both AR and HL (Fig. 8) and testing three different transition densities between snow and firn. Profiles are reported as steps, where each step corresponds to a firn layer. Focusing on the CG15 ice core, we modelled, with both versions, lower densities in the first 4–5 m, with a layer with a particularly low density not matched by the ice core data at around 1–2 m from the surface. This marked decrease however is reduced when a ρ_D of 410 kg m⁻³ is chosen. For higher depths, the model with AR implemented underestimates CG15 density, while with HL a better match of the profile is observed. The best performance is obtained implementing HL and selecting ${\mathit{\rho }}_{\mathrm{D}}=\mathrm{410}\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$, with a mismatch only in the first metres.

Moving to the KCC ice core, the model with AR implemented results in an overestimation of density up to a depth of around 4 m and an underestimation below it. Implementing HL, the density is instead overestimated except for a layer at around 8 m of depth.

Figure 8 Observed and modelled (with multi-layer version) firn density profiles for KCC (b, d) and CG15 (a, c). AR and HL stand for the model versions with Arnaud et al. (2000) and Herron and Langway (1980) implemented, respectively.

4.5 Comparison between multi- and single-layer model versions

In order to understand the approximation introduced modelling firn as a single layer instead of a multi-layer column, we compare in Figs. 9–10 the average firn density obtained with the single-layer version of the model or averaging the density of each individual firn layer obtained with the multi-layer version, weighted for their heights. The results for ${\mathit{\rho }}_{\mathrm{D}}=\mathrm{300}\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$ are not reported since they were not significantly different from the ones with ${\mathit{\rho }}_{\mathrm{D}}=\mathrm{360}\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$. Setting ρ_D to 360 kg m⁻³, we have, implementing HL, a maximum difference between the two average densities of 16.7 kg m⁻³, obtained for the KCC ice core, and, implementing AR, a maximum difference of 7 kg m⁻³, obtained for the CG15 ice core. Higher differences are obtained moving to ${\mathit{\rho }}_{\mathrm{D}}=\mathrm{410}\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$, with a maximum difference of 29 and 14 kg m⁻³ implementing HL and AR, respectively, obtained for the KCC ice core. In all cases the multi-layer version predicts higher average densities.

Figure 9 Average firn density of CG15 ice core observed and obtained with the single-layer version of the model or averaging the density of each individual firn layer obtained with the multi-layer version, weighted for their heights. Each point corresponds to a run of the model whose starting date is reported on the x axis. The ending year for all runs is 30 September 2015. AR and HL stand for the model versions with Arnaud et al. (2000) and Herron and Langway (1980) implemented, respectively. The values 360 and 410 stand for the chosen value of ρ_D.

Figure 10 Average firn density of KCC ice core observed and obtained with the single-layer version of the model or averaging the density of each individual firn layer obtained with the multi-layer version, weighted for their heights. Each point corresponds to a run of the model whose starting date is reported on the x axis. The ending year for all runs is 30 August 2013. AR and HL stand for the model versions with Arnaud et al. (2000) and Herron and Langway (1980) implemented, respectively. The values 360 and 410 stand for the chosen value of ρ_D.

5 Discussion

5.1 Steady-state firn densification

Figure 4 shows a variable performance of the firn densification model depending on the ice core considered; with the exception of KCC and CG15, which show a very good and a very poor performance, respectively, for all the other cores we have a good fit up to a density of about 600–700 kg m⁻³. Bréant et al. (2017), who modified the original model of AR, also observed a variable agreement between the data and model, also for sites with similar accumulation and temperature. They suggested that this may be due to different flow regimes of the sites, since their 1D model does not include this effect. Another consideration that emerges from Fig. 4, also pointed out by Bréant et al. (2017), is that the modelled profile results in worse performances when the observed density profile does not show a clear change in the densification rate near the critical density D₀. The transition is, in fact, more evident for the KCC ice core, which is associated with the best fit. Finally, Bréant et al. (2017) reported a tendency of the model to overestimate the densification rate for lower densities and to underestimate it for higher densities. This is coherent with the results obtained, in which HL predicts lower densities before D₀ and higher densities after D₀ if compared with AR. In order to compare modelled and observed profiles in Fig. 4, it is important to point out that the two models assume stationary conditions; therefore they are not able to reproduce possible changes in the glaciological characteristics. In some of the ice cores it is possible to see a bend in the profile in correspondence to about 20–30 m. The reason could be a combination of ice flow and the upstream effect, i.e. changes in snow accumulation upstream, and these effects cannot be reproduced by a 1D model like the ones used.

5.2 Snow accumulation

Snow accumulation at CG is characterized by a high spatial variability (Keck, 2001; Licciulli et al., 2020). The difference in net annual accumulation of CG11 and CG15, which are about 50 m apart, ranges from $+\mathrm{0.13}\times {\mathrm{10}}^{\mathrm{3}}$ to $-\mathrm{0.266}\times {\mathrm{10}}^{\mathrm{3}}$ $\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ in the period 2002–2012, while the one of CG15 and KCC, which are about 120 m apart, ranges between $+\mathrm{0.41}\times {\mathrm{10}}^{\mathrm{3}}$ and $-\mathrm{0.15}\times {\mathrm{10}}^{\mathrm{3}}$ $\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$. Given the high variability in the accumulation rate, three ice cores may not be enough to fully represent the site. In addition, ice core data are biased due to the fact they are drilled preferentially in the north flank, where accumulation is lower. While the modelled average annual accumulation is in the range of the ones estimated for the north flank of CG (Licciulli et al., 2020), the model is not able to reproduce the observed spatial variability. Due to the lack of dependence on topography in the presented model, we do not expect the model to correctly follow one or the other core. The topography influences the amount of solar radiation received, which in turn influences melting and wind erosion. At the same time, the topography modifies wind speed, which in turn modifies topography itself. This results in a quasi-random spatial variation and a systematic temporal variation in surface accumulation at a given location (Keck, 2001). Surface snow temperature was set equal to air temperature, instead of solving the full surface energy balance that would have required a higher availability of data; surface temperatures, in fact, may also reach 0 ^∘C for air temperatures below 0 ^∘C, mainly when calm conditions are present, or, on the contrary, melting may not occur during positive air temperatures, particularly when wind is present (Keck, 2001).

From the box plots in Fig. 6 we observe that the conserved snow is made up mainly of summer precipitation and that the conserved fraction of solid precipitation reflects the number of hours with greater-than-zero temperature rather than the seasonality of precipitation. The accumulation is, in fact, mainly governed by the wind erosion (Wagenbach et al., 1988) and the presence of wet layers or ice crusts, as well as a faster compaction when temperatures are higher, which protects snow from wind erosion. This is well known in ice core studies at CG, since it results in isotope records that are biased towards precipitation of the warm season (Schöner et al., 2002; Bohleber et al., 2013, 2018). To our knowledge, no models have been applied at CG that try to model the influence of wind speed and temperature on snow accumulation. The only confirmation of this behaviour we have, other than from ice core analysis, comes from temporary measurements of the snow height in nearby sites or observations at CG. For example, at Seserjoch, 4300 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$, the snow height was measured between 1998 and 2000 by Suter et al. (2001), and a main accumulation from about April to November, with practically no accumulation in high winter, was observed. Assessing the link between temperature, wind speed and snow accumulation may be important under scenarios of climate warming for glacial sites, like CG, whose behaviour is strongly regulated by wind activity. At CG, our results suggest that higher temperatures may lead to a counterintuitive response in the snow mass balance, as already conjectured by Alean et al. (1983). In fact, an increase in melting due to a warmer climate may be compensated for, or exceeded, by the consequent reduction in snow erosion, therefore resulting in a net snow accumulation increase.

5.3 Firn density

5.3.1 One-layer model version

The modelled densities have an opposite behaviour when compared with CG15 or KCC. In the first case, the average density is underestimated and in the second case overestimated. The multi-layer version of the model allows us to better understand the reasons for this mismatch. An error in the density of an individual layer will, in fact, affect all the modelled average firn densities that contain that individual layer. This is the case for CG15, where the underestimation of the density of the most superficial layers (Fig. 8) influences the average density of the whole firn column. This influence however decreases when more annual layers are averaged. In the case of CG15, this is probably related to the characteristics of the CG15 location since underestimation was also observed when applying the firn densification model alone in a steady state and not inside the presented snow–firn model. On the contrary, the steady-state application of HL and AR resulted in a good match with the KCC ice core (Fig. 4). Hence, it is reasonable to suppose that the mismatch is related to the modelled snow characteristics. The discrepancy may be explained by the different estimated surface accumulations; the model, in fact, results in an average snow accumulation of 0.46×10³ against the 0.22×10³ $\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ accumulation of KCC. Since the firn densification is driven by overburden stress, this may lead to a systematic bias in the results.

5.3.2 Multi-layer model version

Implementing HL inside the model, we obtained in general a better agreement and a higher variability in the density. The latter is probably related to a higher sensitivity of the equations to surface density, which has instead a low influence on AR as also reported by Bréant et al. (2017). Choosing ${\mathit{\rho }}_{\mathrm{D}}=\mathrm{300}\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$ and looking at the KCC ice core, we observe a decrease in the density at around 7–8 m in both cores. This low density value corresponds to the layer deposited in the 2006–2007 hydrological year. This year was characterized by a solid precipitation event that occurred at the end of the hydrological year. As a consequence, the modelled average snow density decreased from 412 kg m⁻³ at the beginning of 25 September 2007 to 260.6 kg m⁻³ at the end of 27 September 2007. Due to the low amount of snow that is not eroded, new solid precipitation events, in fact, may have a high influence on average snow density. The choice of adopting firn densification equations only after a certain density is reached rather than at the end of the water year allows us to better model this circumstance, in which at the end of the water year the layer has characteristics that are more similar to recent snow than to 1-year-old snow. In comparing the observations and model results, it is also important to point out that, due to the differences in firn height between model and ice core, the ages of the modelled and observed steps do not correspond.

The present model does not resolve a full energy balance to compute surface snow temperature, thus not taking into account the effect of the different amount of solar radiation in the different ice core locations, which is very likely responsible for the significantly different behaviours of the two ice cores. The ice core CG15 is closer to the axis saddle and presumably in a location with a higher exposure to solar radiation than KCC since accumulation is doubled. Consequently, melting events may be more frequent in the CG15 ice core. The differences in the model parameters when run for the KCC or CG15 ice core are the values of MAFT, ρ_D and γ^′ when AR is implemented, and this limits the ability of the model to adapt to the different ice cores.

6 Conclusions

In this study we have proposed a local model that combines snow and firn dynamics. It was derived from the mass balance, momentum balance and rheological equations of snow and firn, combined with semi-empirical and empirical approaches proposed in the literature to model the included mass fluxes. It requires for input hourly (or sub-hourly) series of precipitation, temperature and wind speed with which the series of snow, water and ice inside the snowpack and firn height along with dry-snow and firn density are computed. Two versions of the model were proposed: (1) a version (multi-layer) that considers separately each firn layer and (2) a version (single-layer) that models firn and underlying glacier ice as a single layer. The two implementations allow us to cover different purposes. A simpler model is, in fact, more suitable to reproducing firn inside a hydrological model, where the whole depth–density profile is not necessary and where a reduced number of equations may allow an easier integration. On the other hand, a model that resolves density with depth allows us to assess the influence of meteorological variables on snow and firn characteristics. In order to obtain this, we integrated two existing firn densification models (other firn densification models could possibly be implemented) into a wider model, therefore moving the boundary of the model from surface accumulation and density to hourly meteorological series. While this may not be needed to retrieve past climatological information from ice cores, it is required to assess the response of the system to present and future changes in the climate.

Both models' versions require the calibration of three parameters: a, e and γ^′. In addition, the parameters of the firn densification equation chosen may also need calibration if applied to a temperate site. The modularity of the model allows us to easily test different modelizations of the included fluxes. Depending also on the availability of measured data at the application site, less empirical approaches could be adopted for the mass fluxes. Also, depending on the specific application, some processes presented in the model may be neglected with a reduction in the total number of parameters. A high number of parameters are, for example, associated with wind effects that may be neglected in lower-altitude sites but that are likely to be important in high-altitude sites (Haeberli and Alean, 1985). At Colle Gnifetti, for example, we saw that observed surface accumulations and densities cannot be explained if wind effects are not considered. The modelled reduction in annual snow accumulation due to wind erosion is on average 2×10³ $\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ while the increase in end-of-year snow density is on average 65 kg m⁻³. The strong wind erosion also results in a greater correlation between the amount of snow preserved in each month and the number of days with above-zero temperature rather than with the solid precipitation seasonality. This behaviour may be important in understanding the response of the site to a warming climate. All these elements result in a strong spatial and temporal variability in snow density and accumulation at CG that is not captured by the model. This would, in fact, require taking into consideration the influence of topography on wind speed and the effect of solar radiation on surface snow temperature.

The aim of this study was to illustrate the new snow–firn modelling and to present its potentiality through a case study. In order to integrate it into a hydrological model, further steps are required; in particular fluxes among the different columns should be properly considered regarding both run-off and snow transport. In the present modelling, for example, we considered only the snow erosion. At CG the deposition can be neglected due to its characteristics. The site is a saddle with west–east orientation where the predominant wind direction coincides with the saddle orientation. Besides, an ice cliff at the end of it works as a perfect sink for the eroded snow. A distributed modelling approach needs to take into account the presence of sinks and sources of snow transport. This will also improve the representation of wind effects on snow accumulation since, in order to include snow transport in a 1D framework, approximation of the process is unavoidable for the nature of the process itself.

Appendix A: Complete description and derivation of the snow–firn model

A1 Mass balance equations

The mass balance equations of snow (M_S), liquid water in snow (M_W), refrozen meltwater and rain in snow (M_MF), and firn (M_F) are as follows:

$$\begin{array}{}\text{(A1a)}& {\displaystyle \frac{\mathrm{d}{M}_{\mathrm{S}}}{\mathrm{d}t}}=P_{\mathrm{S}}-M-Q-E_{\mathrm{S}},\text{(A1b)}& {\displaystyle \frac{\mathrm{d}{M}_{\mathrm{W}}}{\mathrm{d}t}}=P_{\mathrm{R}}+M+F-O-E_{\mathrm{W}},\text{(A1c)}& {\displaystyle \frac{\mathrm{d}{M}_{\mathrm{MF}}}{\mathrm{d}t}}=-F-E_{\mathrm{MF}},\text{(A1d)}& {\displaystyle \frac{\mathrm{d}{M}_F}{\mathrm{d}t}}=-O_{\mathrm{F}}+P_{\mathrm{F}}+E_{\mathrm{S}}+E_{\mathrm{W}}+E_{\mathrm{MF}}.\end{array}$$

P_S and P_R are the mass of solid and liquid precipitation events, and they are equal to $P_{\mathrm{S}}=s\cdot {\mathit{\rho }}_{\mathrm{NS}}$ and $P_{\mathrm{R}}=r\cdot {\mathit{\rho }}_{\mathrm{W}}$. The fresh-snow density was calculated as proposed by Liston et al. (2007): ${\mathit{\rho }}_{\mathrm{NS}}={\mathit{\rho }}_{{\mathrm{NS}}_{\mathrm{0}}}+{\mathit{\rho }}_{{\mathrm{NS}}_{\mathrm{w}}}$. Following Anderson (1976), ${\mathit{\rho }}_{{\mathrm{NS}}_{\mathrm{0}}}=\mathrm{50}$ kg m⁻³ if the air temperature $T_{\mathrm{A}}<-\mathrm{15}{}^{\circ }\mathrm{C}$ and ${\mathit{\rho }}_{{\mathrm{NS}}_{\mathrm{0}}}=\mathrm{50}+\mathrm{1.7}\cdot (T_{\mathrm{A}}+\mathrm{15}{)}^{\mathrm{1.5}}$ kg m⁻³ otherwise. The second term gives the increase in fresh-snow density due to wind, and it is computed as ${\mathit{\rho }}_{{\mathrm{NS}}_{\mathrm{w}}}=D_{\mathrm{1}}+D_{\mathrm{2}}(\mathrm{1.0}-\exp (-D_{\mathrm{3}}(u_{\mathrm{2}}-\mathrm{5})\left)\right)$, where u₂ is the wind speed at a 2 m height, D₁=25 kg m⁻³, D₂=250.0 kg m⁻³ and D₃=0.2 s m⁻¹. When u₂≤5 m s⁻¹, ${\mathit{\rho }}_{{\mathrm{NS}}_{\mathrm{w}}}=\mathrm{0}$ kg m⁻³.

M is the snowmelt mass flux that was computed with a temperature-index approach (Hock, 2003): $M=(I\cdot a)(T_{\mathrm{A}}-T_{\mathit{\tau }}){\mathit{\rho }}_{\mathrm{S}}$.

F is the melt–freeze mass flux that was modelled with a coupled melt–freeze temperature-index approach: $F=(I^{*}\cdot e\cdot a)(T_{\mathrm{A}}-T_{\mathit{\tau }}){\mathit{\rho }}_{\mathrm{W}}$.

The run-off O was modelled with a matrix flow approach, and it is equal to O=ρ_WαK_W, where $\mathit{\alpha }={\mathit{\alpha }}^{\prime }(\mathrm{5.47}\times {\mathrm{10}}^{\mathrm{5}}{\mathrm{m}}^{-\mathrm{1}}{\mathrm{s}}^{-\mathrm{1}}$) (DeWalle and Rango, 2008) with ${\mathit{\alpha }}^{\prime }=\mathrm{3600}\mathrm{s}{\mathrm{h}}^{-\mathrm{1}}$. Following Colbeck (1972), K_W was computed as $K_{\mathrm{W}}=K{S}^{*\mathrm{3}}$, in which K is the intrinsic permeability of snow in square metres and S^* is the effective saturation degree of the mixture equal to $S^{*}=(S_{\mathrm{r}}-S_{{\mathrm{r}}_i})/(\mathrm{1}-S_{{\mathrm{r}}_i})$ with $S_{{\mathrm{r}}_i}$ being the irreducible saturation degree equal to $S_{{\mathrm{r}}_i}=\mathrm{0.02}{\mathit{\rho }}_{\mathrm{S}}/\left({\mathit{\rho }}_{\mathrm{W}}\mathit{\varphi }\right)$ (Kelleners et al., 2009) and S_r the average saturation degree of the porous matrix equal to 1 when h_W≥ϕh_S and equal to $h_{\mathrm{W}}/\left(\mathit{\varphi }{h}_{\mathrm{S}}\right)$ otherwise (Avanzi et al., 2015). The intrinsic permeability is obtained using the parametrization proposed by Calonne et al. (2012), and it is equal to $K=\mathrm{3}{R}^{\mathrm{2}}\exp (-\mathrm{0.013}{\mathit{\rho }}_{\mathrm{S}})$, in which R is the equivalent sphere radius. The radius R is defined as $R=\mathrm{3}/\left(\mathrm{SSA}{\mathit{\rho }}_i\right)$, where SSA is the specific surface area in m² kg⁻¹ that was computed by Avanzi et al. (2015) adapting the formula proposed by Domine et al. (2007) with $\mathrm{SSA}=-\mathrm{30.82}\ln ({\mathit{\rho }}_{\mathrm{S}}\times {\mathrm{10}}^{-\mathrm{3}})-\mathrm{20.60}$. When S_r>0.5, to avoid numerical instability, the run-off is calculated with a kinematic wave approximation (De Michele et al., 2013) as ${\mathit{\rho }}_{\mathrm{W}}{\mathit{\theta }}_{\mathrm{W}}{h}_{\mathrm{W}}^{\mathrm{1.25}}$.

The firn melting O_F, which may occur only when the snowpack is absent, was modelled with a temperature-index approach, and it is equal to $O_{\mathrm{F}}=(I_{\mathrm{F}}\cdot a)(T_{\mathrm{A}}-T_{\mathit{\tau }}){\mathit{\rho }}_{\mathrm{F}}\mathit{\delta }\left(h_{\mathrm{S}}\right)$.

P_F is the effect of rain on firn that, when the snowpack is absent, causes an increase in M_F when T_A<0 ^∘C because rainfall is chilled to the firn temperature and a decrease when T_A>0 ^∘C because the energy supplied by rain will be used to melt ice. In the first case $P_{\mathrm{F}}={\mathit{\rho }}_{\mathrm{W}}\cdot r$, while in the second case P_F was set to zero due to its small contribution to mass balance (Doyle et al., 2015).

The terms E_S, E_W and E_MF move the mass of the snowpack still on the ground at the end of each melt season inside the firn, and they are equal to $E_j={\sum }_i{\mathit{\rho }}_j\frac{h_j}{\mathrm{d}t}\mathit{\delta }(t-t_i)$ with j denoting S, MF and W.

Q is the mass of snow eroded by wind obtained from snow transport. The latter was computed adopting the parametrization proposed by Pomeroy et al. (1993) as the sum of a transport in saltation and a transport in suspension.

The saltation transport rate Q_salt ($\mathrm{kg}{\mathrm{m}}^{-\mathrm{1}}{\mathrm{s}}^{-\mathrm{1}}$) occurs only when wind exceeds a given threshold, and it is computed as follows:

$$\begin{array}{}\text{(A2)}& Q_{\mathrm{salt}}={\displaystyle \frac{\mathrm{0.68}{\mathit{\rho }}_{\mathrm{a}}{u}_{\mathrm{t}}^{*}}{u^{*}g}}(u^{{*}^{\mathrm{2}}}-u_{\mathrm{t}}^{{*}^{\mathrm{2}}}),\end{array}$$

where ρ_a is the atmospheric density (kg m⁻³) and u^* and $u_{\mathrm{t}}^{*}$ are the atmospheric friction velocity and the friction velocity applied to the snow surface at the transport threshold, respectively (m s⁻¹). To move from the measured wind speed u to u^*, knowledge of the aerodynamic roughness height z₀ is required. This passage is not straightforward since the value of z₀ during blowing snow events is different from the one during non-transport conditions and it depends on friction velocity (Pomeroy and Gray, 1990). In order to avoid an iterative procedure, we adopted the approximation proposed by Pomeroy and Gray (1990), using $u^{*}\approx \mathrm{0.02264}{u}^{\mathrm{1.295}}$ and $z_{\mathrm{0}}=\frac{\mathrm{0.1203}{u}^{{*}^{\mathrm{2}}}}{\mathrm{2}g}$, where u is the 10 m wind speed (m s⁻¹).

Suspension transport, which occurs only when particles are already in saltation, was computed as follows:

$$\begin{array}{}\text{(A3)}& Q_{\mathrm{susp}}={\displaystyle \frac{u^{*}}{\mathit{\kappa }}}\underset{h^{*}}{\overset{z_{\mathrm{b}}}{\int }}\mathit{\eta }\left(z\right)\ln \left({\displaystyle \frac{z}{z_{\mathrm{0}}}}\right)\mathrm{d}z,\end{array}$$

where Q_susp is in $\mathrm{kg}{\mathrm{m}}^{-\mathrm{1}}{\mathrm{s}}^{-\mathrm{1}}$, κ is the von Kármán constant equal to 0.4, h^* is the lower boundary for suspension equal to $h^{*}=c_{\mathrm{H}}{u}^{{*}^{\mathrm{1.27}}}$ (Pomeroy and Male, 1992) with $c_{\mathrm{H}}=\mathrm{0.08436}{\mathrm{m}}^{-\mathrm{0.27}}{\mathrm{s}}^{\mathrm{1.27}}$, z_b is the top of the surface boundary layer for suspended snow and η(z) is the mass concentration of suspended snow (kg m⁻³) at height z. The mass concentration can be approximate as $\mathit{\eta }\left(z\right)=\mathit{\eta }\left(z_{\mathrm{r}}\right)\exp (-A_Q((B_Q{u}^{*}{)}^{-\mathrm{0.544}}-z^{-\mathrm{0.544}}))$ (Pomeroy and Male, 1992), where η(z_r) is the reference mass concentration for suspension set to 0.8 kg m⁻³ (Pomeroy and Male, 1992), A_Q is equal to 1.55 m^0.544 and B_Q is equal to 0.05628 s^−0.544. z_b was set to 5 m since its value is typically between 5 and 10 m (Déry and Taylor, 1996). The exact value is unimportant because of small mass fluxes at this height (Pomeroy et al., 1993). The snow erosion in the control volume of the model was set equal to the sum of these two transports.

The critical threshold above which snow transport occurs was computed adopting the formula proposed by He and Ohara (2017). The critical shear stress for snow movement can be therefore computed as follows:

$$\begin{array}{}\text{(A4)}& \begin{array}{rl}& {\mathit{\tau }}_t=\\ & {\displaystyle \frac{\begin{array}{c}\left(\mathrm{8}R\cdot C_{\mathrm{g}}\cdot g\right)\left({\mathit{\rho }}_{\mathrm{S}}-{\mathit{\rho }}_{\mathrm{a}}\right)\cos \left(\mathit{\pi }/\mathrm{3}-S\right)+\left(\mathit{\pi }{C}_{\mathrm{c}}\mathit{\varsigma }\right){\left(\frac{C}{R^m}{t}_{\mathrm{d}}\right)}^{\mathrm{2}/n}\\ \times \left(\sin \left(\mathit{\pi }/\mathrm{3}+S\right)+{\left(\frac{C}{R^m}{t}_{\mathrm{d}}\right)}^{\mathrm{1}/n}\right)\end{array}}{\mathrm{2}\left(C_{\mathrm{d}}\sin \left(\mathit{\pi }/\mathrm{3}-S\right)+C_{\mathrm{l}}\cos \left(\mathit{\pi }/\mathrm{3}-S\right)\right)}},\end{array}\end{array}$$

where R is the grain radius (m); t_d is the time since deposition in seconds; C_c, C_d, C_g and C_l are dimensionless coefficients set to 1, 4, $\mathrm{1.3}\mathit{\pi }/\mathrm{6}$ and 3.4 (He and Ohara, 2017); ς is the stress caused by cohesion of ice computed as $\mathit{\varsigma }=\mathrm{1.51}\exp \left(\mathrm{0.44}\right(T_{\mathrm{A}}+\mathrm{9}\left)\right)+\mathrm{6.8}$ (Hosler et al., 1957) for temperatures between −20 and 0 ^∘C with ς in newtons per square metre and T_A in degrees Celsius; and $S=\arcsin \left({\left(\frac{C}{R^m}{t}_{\mathrm{d}}\right)}^{\mathrm{1}/n}\right)$. C, m and n are parameters that influence the rate of ice sintering, modelled following Maeno and Arakawa (2004). In particular, $C=C_{\mathrm{0}}\exp \left(\frac{-Q_s}{R_{\mathrm{G}}(T_{\mathrm{A}}+\mathrm{273.15})}\right)$, in which R_G is the gas constant and T_A is computed as the average air temperature since deposition, $C_{\mathrm{0}}=\mathrm{4.14}\times {\mathrm{10}}^{\mathrm{19}}$ m³ s⁻¹, and $Q_s=\mathrm{1.965}\times {\mathrm{10}}^{\mathrm{5}}$ J mol⁻¹. Finally, m and n are empirical parameters set to 2.9 and 5, respectively, following the results of He and Ohara (2017). Once the critical shear stress is obtained, it is possible to move to the critical friction velocity as follows: $u_t^{*}=\sqrt{{\mathit{\tau }}_t/{\mathit{\rho }}_{\mathrm{a}}}$. If wind speed is lower than the critical threshold, no erosion occurs and Q is set to zero.

To implement the snow erosion routine, we proceeded as explained in the following. When the first solid precipitation event occurs in a time step, the amount of new snow on the ground at the end of the time step, S_A ($\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$), is saved along with the time of deposition and ρ_NS of the event. During the following step, four different situations are possible: (1) a new snow event occurs in the time step. In this case S_A is moved into a vector S_R​​​​​​​ with its time of deposition and ρ_NS. S_A is recomputed as ρ_NS⋅s. (2) T_A>0 ^∘C. In this case S_A and Q are set to 0 $\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$, and all the old snow events memorized in S_R are removed. (3) T_A<0 ^∘C and Q<S_A. In this case S_A is set to $S_{\mathrm{A}}=S_{\mathrm{A}}-Q$. (4) T_A<0 ^∘C and Q>S_A. In this case, if S_R has no elements, Q is set equal to S_A and S_A to 0 $\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$; otherwise the difference between Q and S_A is subtracted from the most recent event in S_R, given that wind speed is higher than the threshold recomputed with the characteristics of that event, and this event is removed from S_R. This is repeated until an event in S_R that cannot be eroded by wind is encountered or the total amount of snow eroded in that time step reaches Q. In the latter case the actual transport is Q, while in the former Q is given by the total amount of snow eroded before reaching the non-erodible layer. The new S_A is the amount of snow associated with the last event considered.

Given that M_j=ρ_jh_j and ρ_k=const with j denoting S, MF, and W and k denoting MF and W, after some algebra we can move from Eqs. (A1a)–(A1d) to Eqs. (2a)–(2d).

A2 Snow densification

The densification of dry snow due to compaction was modelled adopting the formula proposed by Liston et al. (2007):

$$\begin{array}{}\text{(A5)}& {\displaystyle \frac{\mathrm{d}{\mathit{\rho }}_{\mathrm{S}}}{\mathrm{d}t}}=(c\cdot A_{\mathrm{1}}\cdot U){\mathit{\rho }}_{\mathrm{S}}\exp (-B\cdot (T_{\mathit{\tau }}-T_{\mathrm{S}})-A_{\mathrm{2}}\cdot {\mathit{\rho }}_{\mathrm{S}}),\end{array}$$

where $c=\mathrm{0.10}\cdot \mathrm{3600}\mathrm{s}{\mathrm{h}}^{-\mathrm{1}}$, A₁=0.0013 m⁻¹, $A_{\mathrm{2}}=\mathrm{0.021}{\mathrm{m}}^{\mathrm{3}}{\mathrm{kg}}^{-\mathrm{1}}$, $B=\mathrm{0.08}{\mathrm{K}}^{-\mathrm{1}}$ and U is the wind speed contribution (m s⁻¹). For wind speeds $\ge \mathrm{5}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$, $U=E_{\mathrm{1}}+E_{\mathrm{2}}(\mathrm{1.0}-\exp (-E_{\mathrm{3}}(u_{\mathrm{2}}-\mathrm{5.0})\left)\right)$, with E₁, E₂ and E₃ equal to 5.0 m s⁻¹, 15.0 and 0.2 s m⁻¹, respectively, and u₂ is the wind speed at a 2 m height. For wind speed $<\mathrm{5}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$, $U=\mathrm{1}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$. Adding the densification due to mass variation (see De Michele et al., 2013), the total densification rate can be computed as follows:

$$\begin{array}{}\text{(A6)}& \begin{array}{rl}{\displaystyle \frac{\mathrm{d}{\mathit{\rho }}_{\mathrm{S}}}{\mathrm{d}t}}& =(c\cdot A_{\mathrm{1}}\cdot U){\mathit{\rho }}_{\mathrm{S}}\exp (-B\cdot (T_{\mathit{\tau }}-T_{\mathrm{S}})-A_{\mathrm{2}}\cdot {\mathit{\rho }}_{\mathrm{S}})\\ & +{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{NS}}-{\mathit{\rho }}_{\mathrm{S}}}{h_{\mathrm{S}}}}s,\end{array}\end{array}$$

where we assumed that melting events and snow erosion occur at ρ_S=const.

A3 Arnaud et al. (2000) model

The model of AR separates the densification of firn into three stages. The first stage is governed by settling, and it is modelled by Bréant et al. (2017) adapting the equation proposed by Alley (1987). The second stage, which starts when the relative density $D={\mathit{\rho }}_{\mathrm{F}}/{\mathit{\rho }}_i$ equals D₀, is dominated by power law creep, and it is modelled following Arzt (1982) and Arzt et al. (1983). Grains are considered spheres, and each sphere is allowed to increase in radius around fixed centres. Starting from an initial radius l, the new radius l^′ (in units of the initial particle radius l) is $l^{\prime }\left(D\right)=(D/D_{\mathrm{0}}{)}^{\mathrm{1}/\mathrm{3}}$. The growth of spheres increases the number of particle contacts Z from the initial value Z₀ to $Z\left(D\right)=Z_{\mathrm{0}}+b(l^{\prime }-\mathrm{1})$, in which b=15.5. The overlap due to the growth of particles produces an excess volume of material. This excess is distributed uniformly around the portion of the surface of the spheres not in contact. From this excess volume, it is possible to calculate the new radius l^′′ as

$$\begin{array}{}\text{(A7)}& l^{\prime \prime }=l^{\prime }+{\displaystyle \frac{\mathrm{4}{Z}_{\mathrm{0}}(l^{\prime }-\mathrm{1}{)}^{\mathrm{2}}(\mathrm{2}{l}^{\prime }+\mathrm{1})+b(l^{\prime }-\mathrm{1}{)}^{\mathrm{3}}(\mathrm{3}{l}^{\prime }+\mathrm{1})}{\mathrm{12}{l}^{\prime }\left(\mathrm{4}{l}^{\prime }-\mathrm{2}{Z}_{\mathrm{0}}(l^{\prime }-\mathrm{1})-b(l^{\prime }-\mathrm{1}{)}^{\mathrm{2}}\right)}}.\end{array}$$

The average contact area (in units of l²) can be obtained averaging over all existing contacts:

$$\begin{array}{}\text{(A8)}& \begin{array}{rl}{a}_{\mathrm{c}}\left(D\right)& =a_{\mathrm{c}}\left(l^{\prime \prime }\right)={\displaystyle \frac{\mathit{\pi }}{\mathrm{3}Z{l^{\prime }}^{\mathrm{2}}}}\left(\mathrm{3}({l^{\prime \prime }}^{\mathrm{2}}-\mathrm{1})Z_{\mathrm{0}}\right.\\ & \left.+{l^{\prime \prime }}^{\mathrm{2}}b(\mathrm{2}{l}^{\prime \prime }-\mathrm{3})+b\right).\end{array}\end{array}$$

The value of Z₀ for a given value of D₀ is obtained, as proposed by Arnaud et al. (2000), assuming that the effective stress $P^{*}=\left(\mathrm{4}\mathit{\pi }P\right)/\left(a_{\mathrm{c}}ZD\right)$ approaches P as D tends to 1. The third stage begins when pores start becoming isolated (D>D_c) and densification is calculated considering the deformation of ice shells surrounding cylindrical pores (Wilkinson and Ashby, 1975). As for Eq. (2e), the total densification rate is obtained adding the densification due to new mass addition (Eq. 2f).

Appendix B: List of symbols and abbreviations used

Table B1 List of the symbols and abbreviations used in the mass fluxes of the snow–firn model with the exclusion of snow erosion (from A to S).

Table B2 List of the symbols used in the mass fluxes of the snow–firn model with the exclusion of snow erosion (from T to Z and Greek letters).

Table B3 List of the symbols used to compute snow erosion.

Table B4 List of the symbols used in AR.

Table B5 List of the symbols used in HL.

Table B6 List of the symbols related to snow densification rate.

Table B7 List of other symbols and abbreviations used throughout the paper.