Water flow in a natural snow cover is generally a complex process because of the strongly stratified and changing structure of the snowpack. Important differences in, for example, density and grain size between layers cause sometimes very sharp transitions in hydraulic properties. Experiments have shown that water accumulating on capillary barriers in the snow cover can reduce its strength and thereby favour wet snow avalanche formation. Refreezing at capillary barriers within the snow cover leads to the formation of ice lenses or ice crusts, which favors the development of large-grained weak layers (facets and depth hoar). During the melt season, crusts and large grains strongly influence the vertical water flow in snow. Blocking of the vertical flow can lead to lateral flow or to delayed arrival of melt water at the bottom of the snowpack. We extended the model SNOWPACK with a solver for 1D unsaturated flow (Richards Equation), which considerably improved melt water runoff estimations. Here we show that Richards Equation is also capable of reproducing capillary barriers at interfaces between snow layers of differing properties; a real modelling challenge. Moreover, the method to determine the hydraulic conductivity at the interface nodes between snow layers plays a major role in the behaviour of simulated water flow. Comparing the arithmetic with the so-called Darcian averaging approach for internodal hydraulic conductivity, we show that the latter is able to reproduce important accumulation of liquid water at interfaces within the snowpack. Refreezing of this water leads to layers that can be considered crusts or ice lenses, with higher densities than layers above and below. We compare the simulation results with observed snow profiles. The ability of SNOWPACK to reproduce melt-freeze crusts and ice lenses will be helpful in future studies on travel times of water through the snowpack or assessing snowpack stability.