Abstract. The physical processes of heat conduction, liquid water percolation, and phase changes govern the transfer of mass and energy in snow. They are therefore at the heart of any physics-based snowpack model. In the last decade, the use of Richards' equation has been proposed to better represent liquid water percolation in snow. While this approach allows the explicit representation of capillary effects, it can also be problematic as it usually presents a large increase in numerical complexity and cost. This notably arises from the problem of applying a water retention curve in a fully-dry medium such as snow, leading to a divergence and degeneracy in Richards' equation. Moreover, the difficulty of representing both dry and wet snow in a single framework hinders the concomitant solving of heat conduction, phase changes, and liquid percolation. Rather, current models employ a sequential approach, which can be subject to non-physical overshoots. To treat these problems, we propose the use of a regularized water retention curve (WRC), that can be applied to dry snow. Combined with a variable switch technique, this opens the possibility of solving the energy and mass budgets in a fully consistent and tightly coupled manner, whether the snowpack contains dry and/or wet regions. To assess the behavior of the proposed scheme, we compare it to other implementations based on loose-coupling between processes and on the state-of-the-art strategies in snowpack models. Results show that the use of a regularized WRC with a variable switch greatly improve the robustness of the numerical implementation, consistently allowing timesteps greater or equal to 900 s, which results in faster and cheaper simulations. Furthermore, the possibility of solving the physical process in a fully-coupled and concomitant manner results in a slightly reduced error level compared to implementations based on the traditional sequential treatment. However, we did not observe any numerical oscillations and/or divergence sometimes associated with a sequential treatment. This indicates that a sequential treatment remains a potentially interesting tradeoff, favoring computational cost for a small decrease in precision.