The Galerkin finite element method was applied to solve the equations governing the physical processes of snowmelting and infiltration of meltwater in a subfreezing snowpack. The numerical model formulation is presented. The solution domain was discretized into 3-node triangular elements. However, the snowpack melting was calculated by 2-node linear elements accompanied by the corresponding 4-node quadrilateral elements along the melting surface within the same meshing system. Equations which govern the surface melting, the water flow, the local and the global heat conductions were coupled and solved within the same time step. The nodes on the moving snow surface were allowed to move to allow for the snowpack thinning during melting. The hydraulic and thermal properties of snow were updated based on the snow bulk density change due to meltwater refreezing. Various simulations were presented and compared for the purpose of identifying dominant effects of the physical system under complex snowmelting scenarios. The model is able to simulate the growth of thin ice crusts near the snowpack surface. The sensitivity analysis shows that the cold content of snow is the major influencing factor which controls the space-time behavior of water infiltration in a subfreezing snowpack.