The transport of meltwater through porous snow is a fundamental process in hydrology that remains poorly understood but essential for more robust prediction of how the cryosphere will respond under climate change. Here we propose a continuum model that resolves the nonlinear coupling of preferential melt flow and the nonequilibrium thermodynamics of ice-melt phase change at the Darcy scale. We assume that the commonly observed unstable melt infiltration is due to the gravity fingering instabililty, and capture it using the modified Richards equation that is extended with a higher-order term in saturation. Our model accounts for changes in porosity and the thermal budget of the snowpack caused by melt refreezing at the continuum scale, based on a mechanistic estimate of the ice-water phase change kinetics formulated at the pore scale. We validate the model in 1D against field data and laboratory experiments of infiltration in snow and find generally good agreement. Compared to existing theory of stable melt infiltration, our 2D simulation results show that preferential infiltration delivers melt faster to deeper depths, and as a result, changes in porosity and temperature can occur at deeper parts of the snow. The simulations also capture the formation of vertical low porosity annulus known as ice pipes, which have been observed in the field but lack mechanistic understanding to date. Our results demonstrate how melt refreezing and unstable infiltration reshape the porosity structure of snow and impacts thermal and mass transport in highly nonlinear ways, which are not captured by simpler models.

Water stored in snow and ice counts for almost 70% of Earth's freshwater volume (Gleick, 1996). Reliable predictions of the hydrological cycle in cold environments, such as terrestrial snowpack and glaciers, remain challenging but are necessary to improve both water resources and geohazards management under climate variability. A fundamental process that remains poorly understood is how surface-generated melt-water released from its frozen state due to heating of the snow-transports and distributes within the snowpack before entering the groundwater or surface water systems. A robust model for meltwater flow through snow is crucial to formulate reliable predictions in larger-scale models of snow cryohydrology and glaciology.One key challenge of modeling snowmelt hydrology is the ability to robustly capture the infiltration rate and storage location of meltwater within the snowpack. While the generation of melt at snow surface can be relatively uniform in space, meltwater infiltration through the underlying snowpack is known to be highly heterogeneous in nature, forming (a) vertical preferential flow pathways that channelize meltwater (e.g., ice pipes) and (b) lateral flow pathways guided by horizontal low permeability zones (e.g., capillary barriers or ice lenses). Both types of preferential pathways have been observed in the field directly or indirectly (