The numerical simulation of water infiltration in heterogeneous porous media is frequently hindered by strong nonlinearities and mathematical degeneracies inherent to the Richards equation, particularly in the fully saturated and nearly dry regimes. Traditional pressure-based or saturation-based formulations often require variable switching or empirical regularization to handle these limits, which can introduce numerical instabilities and mass balance errors. In this work, we present a robust, structure-preserving finite element framework based on a novel, bounded auxiliary variable, $u$, defined through the incomplete Beta function. This transformation ensures that all terms in the weak formulation remain bounded across the entire saturation range, $S \in [0, 1]$, effectively acting as a physics-based preconditioning that eliminates the need for arbitrary regularization. To address the solver challenges posed by layered and heterogeneous soils, we propose a non-overlapping Schwarz domain decomposition method. The coupled problem is enforced through Robin-type transmission conditions that ensure the continuity of water flux and pressure head across interfaces, enabling the use of efficient parallel sparse direct solvers. We demonstrate the performance of the proposed method using Newton iterations on benchmarks with soil inclusions and layered configurations. Numerical results indicate that the $u$-formulation achieves superior stability, requiring significantly fewer iterations compared to classical regularized schemes. This establishes a scalable and robust approach for high-fidelity modeling of nonlinear flow in complex environments.