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.