Variably saturated flow problems in porous media are of significant interest in hydrogeology. This is due not only to their wide range of applications but also to the numerical challenges posed by strong nonlinearities. Accurately resolving saturation fronts within these flows remains a challenge, particularly in cases where strong moisture gradients develop. Various numerical techniques have been proposed to address this phenomenon; however, some methods suffer from instabilities when steep saturation transitions occur, especially when the mesh resolution is insufficient, leading to unphysical oscillations in numerical outputs. In this presentation, we discuss numerical aspects of an adaptive stabilized approach for solving Richards' equation, focusing on accurately capturing saturation fronts while reducing numerical instabilities. Our method dynamically resolves saturation gradients using an adaptive mesh refinement strategy. The framework is based on an adaptive stabilized formulation (\cite{Calo2020, Giraldo2023,Giraldo2023a}) that incorporates a solution-dependent permeability coefficient, allowing the method to refine the mesh in regions where steep moisture fronts develop. This approach not only enhances numerical stability but also provides an on-the-fly error estimator, guiding refinement in regions where instabilities may arise.