We present onePhaseMixtureFoam, an OpenFOAM-based solver for unsaturated flow in porous media, designed to overcome limitations of classical Richards-type models [1]. The solver is based on a hyperbolic momentum formulation derived from mixture theory and adopts a single-phase Navier-Stokes approach that incorporates inertial and relaxation effects. This enables the model to capture non-equilibrium wetting dynamics such as saturation overshoot, sharp infiltration fronts, and gravity fingering, without relying on dynamic capillarity terms or artificial dispersion. The governing equations are discretized via the finite volume method with a modified PIMPLE algorithm. Constitutive relationships for capillarity and permeability follow van Genuchten–Mualem models. We demonstrate the solver’s ability to solve 1D, 2D, and 3D problems. Furthermore, the hyperbolic formulation is coupled with a thermal module that describes energy conservation including latent heat effects. As a test case, we simulate water infiltration into subfreezing snow, a nonlinear problem involving preferential flow, local refreezing, and evolving permeability. We assess the ability of the simulations to reproduce key features such as flow preferential paths, porosity reorganization, and secondary fingering, aligning with experimental and theoretical findings [2]. Overall, this work highlights the advantages of hyperbolic formulations for modeling instabilities in multiphase porous flows. This solver is fully open-source and allows for the implementation of additional physics, such as porosity evolution, thermal non-equilibrium, or reactive transport, while providing a robust platform for high-fidelity simulations of infiltration processes in geophysical and environmental systems.