Non-equilibrium models are applicable in various physical situations, including phase transitions, hysteresis, and chemical reactions, among others. To model the dynamics of such phenomena, partial differential equations are employed with source terms, as seen in Eq. (1.1). This work focuses on situations where we connect states in equilibrium while allowing for non-equilibrium times. In our model, we model non-equilibrium through relaxation. The source term, in this case, acts as an attractor, and the model’s dynamics allow only a small-time far from equilibrium. We use this theory to model a two-equation system in a two-phase fluid in porous media, where different permeabilities of the flow lead to hysteresis phenomena. We show several examples and study the Riemann solution to gain a better understanding of the solution. We investigate the interplay between nonlinear wave groups of water–oil saturation with different permeabilities in the non-wetting phase and the high-contrast heterogeneity of the porous medium in two-space dimensions. This phenomenon is described by a coupled set of time-dependent partial differential equations of hyperbolic–parabolic–elliptic mixed type. We present and discuss a new projection method with relaxation to obtain 1D-analytical Riemann solutions in an idealized homogenized medium, and we use this solution to validate our numerical method. The primary components of obtaining these solutions are shocks, rarefactions, and bifurcations. We obtain 2D-numerical solutions using an operator splitting approach based on physics, highlighting the discretization methods used. We provide robust analytical–numerical examples to verify the theory and illustrate the approach’s capabilities. In particular, we discuss the 1D-2D behavior of wave groups in situations driven by water–oil infiltration problems of permeable reservoir rocks.