Wet‐Snow Metamorphism Drives the Transition From Preferential to Matrix Flow in Snow

In order to explain the poorly understood transition between preferential and matrix flow in snow, we compared observations from a cold‐laboratory experiment with predictions from a multidimensional water transport snow model. We found a good agreement between the modeled and observed evolution of grain size distributions if two or three dimensions are considered by the model, which validates existing theories of snow grain growth. Furthermore, the model reproduced the spatial migration of preferential flow paths with time and a progressive homogenization of snow wetness and structure only if grain growth was simulated. Spatially varying grain growth thus drives the transition from a preferential flow to a matrix flow regime. This transition is faster when grain size and density are lower, or infiltration rates are higher. This explains why preferential flow is more persistent in firn than in snow.

By coupling wetting front propagation with grain growth , wet-snow metamorphism further complicates the understanding of water flow through snow, because snow structure is not invariant with time. Even though models for wet snow metamorphism have been developed based of experiments on real snow samples , these models have rarely been validated for natural snowpacks (Hirashima et al., 2010), and their accuracy and impact on liquid water flow processes through snow remain largely unknown. Here, we first validate wet-snow metamorphism patterns derived using a widely used parametrization of wet-snow metamorphism  against observations from a cold-laboratory experiment that was representative of seasonal snow . We then use the comparison between model simulations and observations to investigate the role that grain-growth rates play in ruling the evolution of water flow regimes with time and in particular the transition from preferential flow to matrix flow.

Methods
The cold-laboratory experiment considered in this paper includes micro-CT observations of wet-snow metamorphism in a nature-identical snow sample (Schleef et al., 2014) that was subjected to periodical melt cycles with no refreezing (see Figure S1 in the supporting information for the schedule of melt cycles and the estimated height of the snow sample throughout the experiment). More details about this experiment are reported in Avanzi et al. (2017) as Experiment M2. We replicated this experiment using a multidimensional water transport model (Hirashima et al., , 2017. This model simulates water flow through snow using the Darcy-Buckingham law and snow-specific parameterizations of the water retention curve  and permeability (Calonne et al., 2012).
We considered three simulation geometries following the most frequent approaches in this field: 3-D (the original geometry of the model by Hirashima et al., 2014), 2-D (similar to the model by Leroux & Pomeroy, 2017), and 1-D (similar to the geometry of the SNOWPACK model, see Wever et al., 2014). In the 3-D version of the model, the simulated volume was set to 10 × 10 × 25 cm (length × width × height) as a compromise between replicating the real geometry of the experiment (50 × 50 × 30 cm) and reducing computational efforts. In the 2-D version, the simulated area was 200 × 25 cm (width × height) to ensure the same number of simulation pixels for each horizontal plane as in the 3-D simulation (20 × 20 pixels). In the 1-D version, the simulated height was 25 cm. The voxel-pixel size was set to 0.5 cm in all simulations to ensure numerical convergence and to match the observed size of preferential flow path size in Katsushima et al. (2013).
The initial and boundary conditions of the model include snow density, classical grain size (g s , in mm; see Fierz et al., 2009 for definition), volumetric liquid water content (LWC, in m 3 m −3 ), and the external water supply rate. All these conditions were specified following measurements taken during the cold-laboratory experiment: Initial grain size was set to 0.1744 mm, density was set to 227 kg m −3 , LWC was set to 0 m 3 m −3 , and the external melt rate was set to 0.66 mm hr −1 . The latter was estimated based on the observed decrease in snow depth during the experiment and the corresponding snow density. The model does not explicitly consider snowmelt, which was thus represented by an external water input at the top of a snow cover of constant height. Again based on the experimental protocol, this periodical input rate was supplied for 7 hr during each melt cycle, while no input rate was applied for 17 hr. As in Hirashima et al. (2014), heterogeneity was set to 20% of mean grain size based on experimental results in Katsushima et al. (2013).
The boundary condition at the bottom of the simulated volume was set to a maximum LWC of 0.05 m 3 m −3 . This condition is the result of a trial-and-error procedure and ensures that no ponding at the base of the simulated volume will occur. Since this boundary condition will bias estimates of snowpack runoff at the bottom of the simulated volume, we estimated snowpack runoff as the vertical, downward flux at a height of 5 cm from the base of the simulated volume.
Wet-snow metamorphism was simulated using the commonly used parametrization for wet-snow grain growth proposed by Brun et al. (1989). We performed both a simulation with and one without grain growth to quantify the impact of this process on water flow through snow. In the following, we will mostly focus on simulated results at a middle height (12.5 cm), which is comparable to the height at which most of the samples were collected in the cold-laboratory experiment (see Avanzi et al., 2017).
To simulate preferential flow, the model relies on a parameterization of the water entry suction, a threshold suction value allowing water to enter a dry porous medium (Baker & Hillel, 1990): If the water content in a voxel was smaller than the residual water content, suction did not follow a water retention curve but was assumed equal to water entry suction (see Hirashima et al., 2014). The water entry suction in snow has been first parameterized by Katsushima et al. (2013) and Hirashima et al. (2014) starting from the equation proposed for soils by Baker and Hillel (1990). The experiments used by Katsushima et al. (2013) and Hirashima et al. (2014) were, however, performed using a relatively high density and grain size compared to those by Avanzi et al. (2017) (density ≥ 386 kg m −3 and grain size ≥ 0.230 mm) and may not be valid for the type of snow considered here. Indeed, we obtained very large values of water entry suction by directly using the parameterization by Katsushima et al. (2013) and Hirashima et al. (2014) for the fine grained snow in this study. Consequently, preferential flow did not form, in contrast with experimental observations. We therefore adapted the parametrization used by Katsushima et al. (2013) and Hirashima et al. (2014) by fitting a quadratic approximation to the data of Katsushima et al. (2013) (Figure S2). This approximation allowed preferential flow to form due to saturation overshoot but should be regarded as a tentative method at this stage due to a lack of observations of suction profiles for the fine snow considered here.
The impact of this modified parameterization on LWC predictions was verified by simulating water infiltration into very fine snow with a grain size of 0.23 mm, comparable to sample S S in Katsushima et al. (2013). Figure S3 shows that except for the saturation-overshoot behavior at the initial wetting front with the quadratic parameterization, the steady state liquid water content is identical for both parameterizations. Saturation overshoot has been associated with gravity-driven instability (DiCarlo, 2013), meaning that the model result from applying the quadratic approximation is consistent with the formation of preferential flow in fine snow observed in Avanzi et al. (2017). Note that in the rest of this study, the 1-D version does not include preferential flow nor a water entry suction (in contrast with the simulations shown in Figure S3) and only solves for Richards equation (matrix flow) consistently with classical matrix flow snow water transport models.

Results
Both the model and the measurements showed a progressive increase in grain size during controlled melting, but only the 2-D and the 3-D models correctly captured the full range of grain sizes observed during the experiment (Figures 1a-1c). The 1-D version of the model systematically overestimated average grain growth (especially during the early stages of melt), likely because it simulated a consistently larger average LWC compared to the 2-D and 3-D approaches. Indeed, Figure S4 shows that horizontally averaged LWC in the 1-D simulation exceeded the residual water content (0.024 m 3 m −3 ) at all heights starting from the third day, while the 2-D and 3-D versions of the model never reached this vertically homogeneous conditions. In contrast, the 2-D and 3-D versions of the model show LWC values close to saturation inside preferential flow paths (see profiles of horizontally maximum LWC in Figure S5). These flow paths were, however, surrounded by dry snow (see Figure 2), such that (1) the average grain growth in the 2-D and 3-D model versions was much slower than in the 1-D version and (2) maximum grain size in the 2-D and 3-D model version was slightly but consistently larger than average grain size in the 1-D model. The mean bias and the Root Mean Square Errors (RMSE) for mean grain size were 0.256, 0.049, and −0.092 mm and 0.32, 0.19, and 0.20 mm for the 1-D, 2-D, and 3-D model versions, respectively.
The temporal derivative of grain size (a measure of grain-growth rate) showed different trajectories based on the number of dimensions considered by the model (Figures 1d-1f). In the 1-D simulation, this derivative decreased with time according to an approximately exponential rate (see Figure 1d). In both the 2-D and  3-D simulations, this derivative was approximately linear, with a negative and a positive slope in the 2-D and the 3-D simulation, respectively (−1.45 × 10 −3 mm day −2 and 1.70 × 10 −3 mm day −2 in 2-D and 3-D, respectively). Thus, the 1-D model version implicitly neglecting preferential flow not only overestimated grain growth but also misrepresented its temporal pattern.
LWC distributions in the 3-D simulation (here assumed as the most representative of observations) showed that the number of preferential flow paths increased through time, thus leading to progressive spatial homogenization of snow wetness and a transition from a preferential flow regime to a matrix flow regime ( Figure 2). This transition occurred after ∼10 days, when most of the upper portion of the sample was homogeneously wet, congruent with matrix flow (Figures 2e and 2f). Grain growth emerged as the key driver of this mechanism because grain growth, induced by wet-snow metamorphism, progressively increased pore size, thus reducing suction in the preferential flow paths. Due to the lower suction of coarser wet snow than that of the adjacent, finer dry snow, some portion of dry and fine snow became easier to infiltrate than the-now coarser-previously formed preferential flow paths (see Figure S6). This feedback mechanism thereby leads to spatial migration of preferential flow paths, which ultimately homogenizes both snow structure (because of coarsening) and wetness. A simulation of the same experiment without grain growth showed that in that scenario water tended to always infiltrate into the same path, with no spatial migration or homogenization ( Figure S7). Validating this no-growth scenario in reality would likely require experiments with alternative liquids and under subfreezing conditions, which we suggest for future work.
In addition to water distributions within the snowpack (Figure 2), the spatially varying grain growth and its effect on liquid water flow also impacted snowpack runoff timing and amount (Figure 3). Instantaneous runoff for the 3-D simulations with and without grain growth were similar during the first day of the simulation, even if the peak according to the simulation with grain growth occurred about 40 min earlier and was about 11% larger than the peak runoff in the simulation without grain growth (Figures 3a-3d). Peaks in instantaneous runoff for the simulation with grain growth then slowly decreased with time, while the simulation without grain growth showed a sharp decrease between the first and the second melt cycle, after which the peaks in runoff remained fairly constant during all the following melt cycles. This contrasting behavior can be explained by the interaction of spatially varying grain growth with preferential flow path formation. When preferential flow paths form, liquid water concentrates at front due to saturation overshoot (Katsushima et al., 2013;Leroux & Pomeroy, 2019), which causes a spike in runoff once the water arrives at the bottom of the snow. While this process is repeated in the simulation with grain growth because new paths are created during each melt cycle, the same preferential flow paths are active during all melt cycles in the simulation without grain growth ( Figure S7), and thus, saturation overshoot only occurs during the first wetting. This implies that, after adjusting suction and conductivity to the external input rate, runoff peaks in the simulation without grain growth become consistent with each other. On the other hand, as the grain size in the whole sample increases (Figure 2), the simulation with grain growth tends to become similar to the simulation without grain growth simply because saturation overshoot decreases and water flow is spatially more homogenously distributed. Both conditions are hydraulically similar to the conditions in the persistent single preferential flow path of the simulation without grain growth.
Cumulative runoff (Figure 3e) exhibits different behavior as well whether or not grain growth is considered. In particular, the simulation with grain growth predicted a lower cumulative runoff than the one without grain growth. This may seem counterintuitive as grain growth coarsened snow structure and as such should have made snow more conductive and less retentive than a simulation without grain growth. However, grain growth also promoted redistribution of water within the snow sample (see Figure 2) and thus translated into a less efficient transport mechanism because some water remained within the sample as residual water content rather than being discharged. Interestingly, the 1-D simulation with grain growth underestimated cumulative runoff during the first 6-8 days, but predictions after 2 weeks agreed very well with a 3-D simulation with grain growth (see Figures 3e and S4). We interpret this agreement as again due to the fact that both the 1-D and the 3-D simulation with grain growth ultimately tended toward a hydraulically similar, spatially homogeneous condition in which matrix flow dominated.
Our simulations focused on replicating laboratory experiments and were performed with a relatively low density and small grain size compared to seasonal snow at the onset of the melt season, or to multiyear firn. In order to clarify the applicability of our results to denser snow, we performed a 3-D simulation with initial density and grain size equal to 350 kg m −3 and 0.5 mm, respectively. For consistency reasons, all other initial and boundary conditions as well as the quadratic approximation for water entry suction were the same as the simulation reported in section 2. This simulation showed that the mechanisms observed in both the laboratory experiment and our simulations were still at play in denser snow, including the creation of new preferential flow paths and the deactivation of others (see Figure 4). However, preferential flow path migration was much slower, and the simulation did not show the emergence of matrix flow within the time simulated (15 days). We attribute this to the decrease of grain growth rates with increasing grain size (Colbeck, 1986).
Also, the external melt rate in the original laboratory experiment was relatively small compared to late-spring melt conditions (0.66 mm/hr, which is 4.62 mm/day). To explore the sensitivity of these results to the water input rate, we performed a set of 2-D simulations with a simulated area of 10 × 25 cm, and using larger water input rates. Three-dimensional simulations would have been too computationally intensive for this sensitivity test. We found that with increasing water input rates (0.66, 1.32, 3, 6, and 10 mm/hr), the percentage of wet area at the end of each simulation (15 days) increased to 29.4%, 31.3%, 41.5%, 47.2%, and 51.3%, respectively. The rate of microstructural changes increases with water input rate, which suggests that the time needed for the transition between preferential flow and matrix flow decreases with increasing water input rate.

Discussion
Our 3-D model and micro-CT data showed a very good agreement in terms of snow microstructural evolution (Figure 1), which implies that the existing theory of wet-snow metamorphism (Brun, 1989) provides an adequate simulation of grain growth during periodical snowmelt if a 2-D or 3-D geometry is used. On the other hand, one-dimensional geometries overestimate grain size when applied to conditions that are more representative of natural snow  than the homogeneous conditions of, for example, Colbeck (1986). Since 1-D approaches are still the most popular choice for firn meltwater retention studies (e.g., Marchenko et al., 2017;Reijmer et al., 2012;Steger, Reijmer, van den Broeke, Wever, et al., 2017), wet-snow avalanche assessments (Wever et al., 2018;Wever, Vera Valero, et al., 2016), and snowmelt forecasting (e.g., DeWalle & Rango, 2011;Markstrom et al., 2015;Wever et al., 2017), we suggest that future modeling improvements in this regard should focus on increasing model dimensionality to improve predictions of the structural evolution of wet snow (Verjans et al., 2019).
Coupling between coarsening and flow regimes has been hypothesized for a long time (Schneebeli, 1995), but we demonstrated here that grain growth may also drive the transition from predominantly preferential flow in the early melt phase to predominantly matrix flow in the main melt phase through spatial migration of preferential flow paths and homogenization of snow structure and wetness. The time needed to pass from the first to the second flow regime was relatively short in our simulation (∼10 days), compared to the duration of a snowmelt season. Together with the fact that seasonal snow covers are only up to a few meters deep, the short time needed by this transition could explain why 1-D, matrix flow-based models have been usually successful in simulating snowmelt at seasonal time scale (DeWalle & Rango, 2011;Wever et al., 2014), even if they are inherently biased when it comes to simulating snow-structural evolution in wet conditions (Figure 1). Spatial migration of preferential flow paths also caused a decrease in cumulative runoff compared to a no-growth scenario, which means that the reconstruction of snow structural evolution in wet conditions directly impacts runoff predictions. Field experiments coupling snowpack runoff measurements and internal structural evolution at high spatial resolution would help shed further light on these interactions. These experiments may also help further develop detailed multilayer one-dimensional models like SNOWPACK, which can simulate preferential flow but currently do not consider different microstructural characteristics in the matrix and preferential domains. Based on our results here, we hypothesize that implementing this may likely reproduce at least part of the heterogeneity in grain size found in Avanzi et al. (2017).
Preferential flow has also been observed in firn, for which it is believed to be the predominant water flow pattern (e.g., Humphrey et al., 2012). Preferential flow paths in firn are often well established and of larger size than observed in seasonal snow, such that the term piping has been associated with it (Pfeffer & Humphrey, 1998). The simulation with denser snow (Figure 4) suggests that the reason why piping is so recurring in firn observations at typical infiltration rates may simply be that coarsening in firn is very slow, if not absent. Studies have tried to replicate this process in one-dimensional bucket-type models by reducing the water holding capacity of firn to cause an acceleration in water flow (Marchenko et al., 2017;Reijmer et al., 2012). Our results suggest that such modifications only apply for the cases where they have been developed for.
The multidimensional model used in this paper does not consider hysteresis  and parametrizes water entry suction using an empirical approach developed primarily for water repellent porous media (see section 2). Using water entry suction instead of hysteresis leads, however, to an unrealistic increase of pressure and can potentially cause model instabilities (Leroux & Pomeroy, 2017). Also, Leroux and Pomeroy (2017) pointed out that neglecting capillary hysteresis leads to an overestimation of capillary pressure at the tip of the preferential flow path, where saturation overshoot takes place (DiCarlo, 2013). Furthermore, using a dynamic capillary pressure (Hassanizadeh & Gray, 1993) has been found to simulate saturation overshoot in snow more accurately (Leroux & Pomeroy, 2019) than using a water entry suction. These notions will form the basis of future work, including further validation of our results using newly introduced treatments of water infiltration through snow based on nonequilibrium Richards equation and dynamic capillary pressure (Leroux & Pomeroy, 2019).

Conclusions
We reproduced observations of wet-snow metamorphism that are representative of field conditions using a multidimensional water transport model for snow. The characteristics of grain growth simulated using different geometries were significantly different: In a 2-D and 3-D simulation, the range of simulated grain growth corresponded well with measurements, while average grain growth rate was systematically overestimated using a 1-D simulation. Preferential flow paths in our simulations migrated with time, progressively expanding the area occupied by wet snow and promoting a transition from preferential to matrix flow. A simulation in the same conditions but with no grain growth did not reproduce spatial migration of preferential flow paths, making grain growth the key driver of this transition in flow regimes. We also found that with increasing grain size and density, and consequently lower grain growth rates, preferential flow paths are more persistent in time, which explains why coarse layers of porous ice such as firn on ice sheets generally show more pronounced preferential flow than layers of finer grains as found in seasonal snowpacks. Furthermore, increasing melt rates were found to decrease the transition time from a preferential to a matrix flow regime. The microstructural changes in the snowpack and the associated water flow patterns were also found to modulate snowpack runoff characteristics, such as peak outflow and cumulative runoff, particularly in early stages of snowmelt.