Ice layer formation in the snowpack due to preferential water flow: case study at an alpine site

Ice layers may form in the snowpack due to preferential water flow, with impacts on the snowpack mechanichal, hydrological and thermodynamical properties. This case study at a high-altitude alpine site aims at monitoring their formation and evolution during winter 2017 thanks to the combined use of a comprehensive observation dataset at daily frequency and detailed snowpack modelling. In particular, daily SnowMicroPen penetration resistance profiles enabled to better identify ice layer temporal and spatial heterogeneity when associated with traditional snowpack profiles and measurements, while 5 upward-looking ground penetrating radar measurements enabled to detect the water front and better describe the snowpack wetting when associated with lysimeter runoff measurements. One-dimensional snowpack simulations with the SNOWPACK model successfully represented the formation of some ice layers when using Richards equation and preferential flow domain parameterization, with a newly implemented ice reservoir. Detailed snowpack simulations with snow microstructure representation, associated with high-resolution comprehensive observation dataset were shown relevant for studying and modelling 10 such complex phenomena, despite limitations inherent to 1D modelling.

. Overview of the WFJ research field, with contour lines in orange and labels indicating the location of the measurements used in this study. of 19.6 mm 2 . Every day during the winter season, five to seven SMP measurements were performed along the corridors, with a 15 cm spacing perpendicular to the direction of the corridor, providing indications about the local spatial heterogeneity of potential ice forms. A limited number of gaps in the daily measurements have to be noted, in particular from 25 March to 27 March and 7 and 8 April in the wetting period.

95
For longer-term observations of ice layers, we gathered traditional snowpack profiles performed every two weeks during the winter seasons 1999/2000 to 2015/2016.

Snowpack simulations 2.2.1 Simulation configuration
SNOWPACK simulations were performed for the WFJ study site in winter 2017. They were initialized with a traditional profile 100 recorded on 3 January 2017, to provide a realistic base of the snowpack with a snow depth of 47 cm. The initialization date was chosen early enough to assess the model ability to simulate the microstructure evolution as well as water percolation. The simulations were driven by optimal in situ meteorological measurements ( Fig. 1) of air temperature and humidity (ventilated sensors), near-surface wind, solar and longwave irradiance (WSL Institute for Snow and Avalanche Research SLF, 2015; Wever 4 https://doi.org/10.5194/tc-2020-24 Preprint. Discussion started: 25 February 2020 c Author(s) 2020. CC BY 4.0 License. et al., 2015). The snowfall input is driven by measured snow depth increments (Wever et al., 2015), enabling direct comparisons to measurements and outcomes of Wever et al. (2016).
Three water transport schemes implemented in SNOWPACK were evaluated. First, the bucket approach (BA) is a common method used in snowpack models (e.g. Bartelt and Lehning, 2002;Vionnet et al., 2012), assuming that water is transported to the next downward layer when the liquid water content exceeds the water holding capacity of a given layer (depending on the ice volumetric content of snow; Coléou and Lesaffre, 1998). Second, the Richards equation (RE) was implemented in 110 SNOWPACK by Wever et al. (2014) to account for capillary effects. These effects are modelled taking into account the water retention curve (van Genuchten, 1980) and the hydraulic conductivity of snow (Mualem, 1976;van Genuchten, 1980;Calonne et al., 2012). Third, a dual-domain approach parameterizing preferential flow and using Richards equation (RE/PF) was recently developed. Water exchanges between the matrix domain and the preferential flow domain are determined according to the water entry pressure head in the matrix layers and the saturation in the preferential flow domain: this implementation is described in 115 details in Wever et al. (2016). The two tuning parameters of this scheme were chosen here accordingly to Wever et al. (2016): the threshold in saturation of the preferential flow domain Θ th = 0.08 and the parameter related to the number of flow paths per square meter N = 0. In particular, N = 0 implies no refreezing of the preferential flow water (Wever et al., 2016). Similarly to Wever et al. (2016), SNOWPACK simulations are carried out at high vertical resolution, with a layer merging threshold of 0.25 cm and new snow layer initialization of 0.5 cm. High resolution is necessary to permit the formation of very thin high density 120 layers.

Implementation of ice reservoir
A new parameterization of ice layer formation due to preferential flow was implemented as a complement to the RE/PF scheme.
It is summarized in Fig. 2. In the RE/PF scheme, when the saturation in the preferential flow domain exceeds the threshold Θ th , water flows back to the matrix domain depending on the freezing capacity of the matrix domain; if the threshold is still 125 exceeded then, saturations in both domains are equalized (Wever et al., 2016). Thus, the part of water that freezes in the matrix domain is spread in the whole layer, while ice lenses may only form locally at the base of the flow fingers.
To overcome this issue, we developed an ice reservoir parameterization. The water normally transferred from the preferential flow domain to the matrix domain that freezes instantly is stored in an ice reservoir (step 4 in Fig. 2), instead of being added to the ice volumetric content of the matrix. The ice reservoir is representative of the volumetric content of ice lenses (i.e. 130 spatially heterogeneous ice) in a given layer. The transferred water that does not freeze goes in the matrix domain , i.e. is spread homogeneously (step 5 in Fig. 2).
Furthermore, the saturation threshold in the PF domain (Wever et al., 2016) was chosen as a simple solution to the unability of Richards equation to model the saturation overshoot present in the tip of flow fingers (DiCarlo, 2007). This simple parameterization can lead to inconsistencies due to the vertical discretization of the simulated snowpack. After water has been 135 transferred to the matrix at the layer corresponding to the finger tip (i.e. where the saturation threshold was exceeded), the highest saturation is then reached more likely at the layer above, where no water transfer occurred. Because of that, the water transfer from PF domain to matrix domain may spread over too many layers, instead of being concentrated in the lowest Figure 2. Scheme of the ice reservoir parameterization in SNOWPACK. Blue represents liquid water, cyan represents ice and red arrows represent water or ice transfers to another domain. Steps 1 to 6 are described in the text. layer (i.e. the tip of the flow finger). To overcome this issue, the ice reservoir was cumulated in the lowest layer. When the ice volumetric content of the cumulated ice reservoir added to the ice volumetric content and water volumetric content of the 140 associated matrix layer exceeds the corresponding ice density threshold of 700 kg m −3 in SNOWPACK, there is enough ice to consider it as horizontally homogeneous: the ice content of the cumulated ice reservoir is then transferred to the associated matrix layer (step 6 in Fig. 2). Simulations with the ice reservoir parameterization are called RE/PF/IceR hereafter.  (Richter et al., 2019). This layer was covered by new snow at the end of December, and several small snow storms lead to a maximum snow depth of 205 cm on 10 March, i.e. slightly lower than the average maximum snow depth. The snowpack had entirely melted on 14 June. Overall, this 150 winter was characterized by lower snow depth than long-term averages in the region. This layer is continuously identified as buried surface hoar (as primary or secondary grain) until the end of March, with a grain size substantially larger than the layer above (consisting of rounded grains or faceted crystals), thus constituting a capillary barrier with a classical fine-over-coarse structure (e.g. Avanzi et al., 2016). Ice is observed above this barrier from 28 March, either as a homogeneous layer or as ice lenses mixed with melt forms (red in Fig. 3). Thicknesses between 0.5 cm and 1 cm were reported. The presence of nearby slopes (NNW of the snow profile corridors in Fig. 1) may suggest a water input through lateral flow along the capillary barrier. Simulations (not taking into account lateral flow) may provide complementary insights to determine whether vertical preferential flow was sufficient to form an ice layer. Layer 2 corresponds to surface hoar 160 appearing in mid-February, and forming a capillary barrier once buried. An ice layer is observed at that level in most profiles after 28 March. of water front could only be derived until 8 April due to technical issues. Before 30 March (first dashed line), no snowpack runoff was observed, the water front was always higher than 1 m, and temperatures in the lowest first meter of the snowpack were all below 0 • C. Between 30 March and 9 April (second dashed line), the water front remained high (mostly higher than 1 m), a small amount of snowpack runoff was observed (less than 3 kg m −2 day −1 ). Snow temperatures gradually increased 170 to reach 0 • C by the end of the period in the lowest meter of the snowpack. From 9 April, all temperatures in the lowest meter reach 0 • C and snowpack runoff increases. Although after 9 April no more water front estimates were available from the upGPR 7 https://doi.org/10.5194/tc-2020-24 Preprint. Discussion started: 25 February 2020 c Author(s) 2020. CC BY 4.0 License. data, it reached the lowest value (85 cm) on 8 April. The snowpack runoff is yet low compared to the more significant runoff starting in mid-May, as shown in Fig. 5: the first 20 kg m −2 of total snowpack runoff are reached on 10 April, while the first day with a daily snowpack runoff higher than 10 kg m −2 day −1 is 13 May.

175
These measurements give insights in the timing of water percolation in the snowpack. Before 9 April, the bottom of the snowpack was cold and dry while the water front was still mostly above 100 cm. The low snowpack runoff values were thus likely due to preferential flow paths reaching the ground. After 9 April, the lowest meter of the snowpack reaches an isothermal state at 0 • C: snowpack runoff increases markedly and is likely due to matrix flow reaching the ground.

Capillary barriers and ice layers 180
To study the formation of ice through preferential flow in subfreezing snow, we focus on the lowest meter of the snowpack.
In this part, the manual snow profiles enable to identify three main capillary barriers (Fig. 3): the two layers of buried surface hoar where ice forms at the end of March (defined as layers 1 and 2 previously) and the top of the depth hoar base layer, where 8 https://doi.org/10.5194/tc-2020-24 Preprint. Discussion started: 25 February 2020 c Author(s) 2020. CC BY 4.0 License.  higher water content is observed in April and May. Layers 1 and 2 are marked by grain size heterogeneity with the overlying layers: on 28 February, 1 mm over 2.5 mm for layer 1, 0.5 mm over 2 mm for layer 2.

185
Daily SMP measurements enable to more clearly identify the temporal and spatial variability of ice formation. Figure 6 represents the evolution of penetration resistance from 1 February to 19 April, with a scale from 0 N to 2 N to highlight variations in dry snow. The deep MFcr is visible in the middle of a low resistance depth hoar layer, at approximately 20 cm.
In March, the buried surface hoar of layer 1 is marked by a lower penetration resistance than the surrounding faceted crystals.

195
The penetration resistance measured at these two layers was tracked in the SMP profiles. To identify these layers, all SMP profiles were superimposed on the traditional profile observed at the closest date, compared to each other for a given day and the next day. The value of SMP penetration resistance was then visually picked. Figure 8 shows an example of this procedure for two SMP profiles of 4 April 2017.
The evolution of the penetration resistance of the two tracked layers is plotted in Fig. 9 from 14 March to 19 April. For 200 layer 1, the penetration resistance remains very low (less than 1 N) until 24 March. It corresponds to the observed layer of buried surface hoar. On 29 March, all seven SMP measurements show penetration resistance higher than 10 N indicating the continuous presence of ice. Afterwards, values alternate between low resistance (less than 5 N) and very high resistance (more than 10 N), as visible in Fig. 8 on 4 April. This is consistent with the visual observations reporting a layer in which pure ice and melt forms are both observed. These observations suggest that water ponding at the capillary barrier did not freeze 205 everywhere on the study plot where we performed these measurements. For layer 2, very low penetration resistances are also measured until 24 March, corresponding to the observed buried surface hoar. After 29 March, values increase (mostly higher than 5 N, often more than 10 N), indicating formation of ice. After 9 April, no more high resistances are measured, rather low resistances corresponding to a wet layer of melt forms. The evolution of penetration resistances for layer 2 show more temporal consistency than layer 1, suggesting that the ice layer disappears totally after 9 April.     Fig. 10a). It is associated with some melting close to the surface leading to preferential water flow refreezing at capillary barrier of layer 1 (Fig. 10b). The water transferred for refreezing in the matrix domain is spread homogeneously which forms a crust with a density of approximately 320 kg m −3 .
The transfer spreads vertically due to the issues mentioned in Sect. 2.2.2. This thick melt-freeze crust was not observed in 220 the manual profiles nor the SMP measurements. Other melt-freeze crusts form close to the surface from mid-February. They were observed (Fig. 3) but were thinner than the simulated ones. A little surface melting is simulated on 1 February and leads to preferential flow (Fig. 10). Contrary to later simulated melting in mid-February, it was not observed: the measured snow surface temperature remained slightly under melting point. This simulation error is likely due to excessive surface turbulent fluxes input. New simulations were run without this melt water input, with no effect on the later snowpack structure due to the 225 limited melting amount. Figure 11 shows the grain shape and the liquid water content in the matrix domain from 15 March to 15 April, i.e. the period of transition from dry to ripe snowpack when ice layers formed (Sect. 3.1). No ice layer forms, except at the snowpack basis, which is probably due to a boundary effect at the interface between snow and soil. The thick melt-freeze crust is still present at layer 1, thus no ice layer forms (Fig. 11a). However, a higher water retention is simulated (Fig. 11b). Due to the 230 excessive formation of melt-freeze crusts, the simulated snow microstructure at layer 2 does not reproduce the observed snow microstructure and the capillary barrier, thus no ice layer forms. Matrix flow reaches the ground and the snowpack is entirely isothermal on 31 March (Fig. 11b), i.e. 9 days before the observations (Sect. 3.1.2).
A sensitivity study on the parameters of the dual-domain approach (Θ th and N ) was performed, but could not resolve the issues about ice formation highlighted here. Simulations were also analyzed in terms of snowpack runoff, confirming 235 earlier findings (Wever et al., 2016;Würzer et al., 2017). The BA scheme underestimates the snowpack runoff, the RE scheme overestimates it, and the addition of preferential flow increases this overestimation (not shown). The onset of snowpack runoff is delayed compared to observations for BA and RE, because preferential flow is not simulated, but RE/PF simulations show the onset of snowpack runoff too early.

With ice reservoir parameterization 240
Simulations were also performed with the ice reservoir parameterization (RE/PF/IceR) to assess its ability to improve the simulation of ice layer formation compared to the previous simulations. Figure 12 shows the grain shape and liquid water content in the PF domain for the month of February. Contrary to the RE/PF simulation, no melt-freeze crust forms at layer 1 ( Fig. 12a): the water leaving the PF domain and refreezing is in too low quantity to be considered as representative of the mean state of the snowpack in this layer, it is thus stored in the ice reservoir. The fine-over-coarse grain transition forming a capillary 245 barrier is preserved. Note that liquid water content in the PF domain (Fig. 12b) is almost not modified compared to the RE/PF simulation (Fig. 10b). Liquid water transport is similar, and in particular the vertical spreading of water ponding, but ice in the reservoir is concentrated at the capillary barrier. At the end of February, less melt-freeze crusts are formed than in the RE/PF simulation, even though the ones surrounding layer 2 are also thicker than observed.   i.e. 4 to 5 days earlier than observed (Sect. 3.1.1 and Sect. 3.1.3). At this date, ice is transferred from the ice reservoir to the matrix domain (Fig. 14). The ice layer formed is 43 mm thick, with a dry density of 821 kg m −3 and a significant volumetric liquid water content of θ matrix = 7.8% and θ P F = 1.7% (on 24 March 13 UTC). Similarly to the RE/PF simulation, no ice 255 forms at layer 2 due to the presence of melt-freeze crusts. However, less melt-freeze crusts are simulated in the snowpack, which is more in accordance with the observations. The matrix water flow reaches the ground on 30 March (Fig. 13b), i.e. 10 days before the observations (Sect. 3.1.2). The ice reservoir does not modify the snowpack runoff compared to RE/PF simulations (not shown).
14 https://doi.org/10.5194/tc-2020-24 Preprint. Discussion started: 25 February 2020 c Author(s) 2020. CC BY 4.0 License.  Wever et al. (2016). Ice layers at the snow-soil interface are not taken into account because of changing snowpack base in the observations over the winters (wooden board, gravel) and possible boundary effects in the simulations. Only simulated ice layers are verified against observations, to calculate hits (number of simulated 265 ice layers matching an observation) and false alarms (number of simulated ice layers that do not match any observation). A height difference of 20 cm is used for matching of simulations and observations, similarly to Wever et al. (2016). We assume the formation date of an observed ice layer is comprised between the last snowpack profile without observed ice layer and the first profile where it is indicated. If the simulation date is more than one month away from the observed formation time interval, or if the height difference is more than 20 cm, the event is considered as a false alarm. Results of this multi-year evaluation are 270 summarized in Table 1.
Overall, the addition of the ice reservoir parameterization to the RE/PF scheme enables to form more ice layers, with a higher number of hits (15 against 6, with a one month tolerance) and a lower number of false alarms (1 against 6). The simulated ice formation date is in average 22 days earlier than the observation interval for the RE/PF scheme and 4 days earlier for the RE/PF/IceR scheme. The RE/PF/IceR configuration also mitigates the number of unobserved early melt-freeze crusts 275 compared to the RE/PF configuration (not shown), as already highlighted for season 2016/2017 ( Fig. 10 and Fig. 12). However, a very high number of ground ice layers were simulated (17 for the RE/PF/IceR scheme against 8 for the RE/PF scheme). This number is probably excessive, and due to possible snow-soil boundary effects.

Discussion
This case study of ice layer formation at Weissfluhjoch enables to assess both a comprehensive observation dataset and detailed 280 1D snowpack simulations for monitoring a complex process. We discuss hereafter the relevance of these methods and results.
First, the combined use of traditional snow profiles with measurements of higher temporal resolution like the SMP provides a suitable observation framework to study the transition period from dry to isothermal snowpack when ice formations appear due to preferential flow. Snowpack runoff measurements associated with snow temperature sensors and upGPR-derived water front gave insights about the homogeneous wetting of the snowpack and the period when the bottom of the snowpack was primarily 285 affected by preferential flow. SMP profiles of penetration resistance showed clear signals of ice presence, when compared to visual observations, with values higher than 10 N (Fig. 7) while penetration resistances in dry snow were mostly lower than 2 N ( Fig. 6) and melt-freeze crusts were characterized by intermediate penetration resistances (usually between 5 N and 10 N, Fig. 7). Identification of ice layers with several SMP profiles regularly spaced also offers a more quantitative estimate of ice heterogeneity than a subjective visual observation. However, the temporal and spatial variabilities may be complex to distin-290 guish, as shown for layer 1 (Fig. 9). The visual layer tracking of SMP profiles is also a source of uncertainties. Hagenmuller and Pilloix (2016) developed a method matching several hardness profiles to synthesize them into one representative profile. This method was not considered relevant for the present study which focuses on local heterogeneity of ice layers. Moreover, when ice layers are too thick, the SMP cannot go through them as happens often in spring. For winter 2017 at WFJ, no complete SMP profile could be performed after 19 April. Finally, the difficulty to identify the exact date of ice formation with biweekly One-dimensional SNOWPACK simulations provide complementary insights to the observation dataset, despite the spatial 300 heterogeneity of ice layer formation due to preferential flow. The addition of an ice reservoir enables to parameterize the local formation of ice at capillary barriers: it may thus be considered as representative of the volumetric content of ice lenses at a