In recent years, there has been a resurgence in developing new solver technologies for addressing highly complex and large-scale flow simulations on specialized parallel and multicore architectures in a very effective manner. Methods such as algebraic multigrid (AMG), Krylov recycling (e.g., deflation, Krylov-secant) and extensions to two-stage preconditioners (e.g., GPR) have been copping the scene in latest solver advances. Nevertheless, there is still a long way to transit to be able to reduce the gap between achieving maximum robustness and parallel efficiency of these solvers in a wide range of problems that the oil industry is currently pursuing on. A new generation of solvers seems to require capabilities to recapture part of the masked physics that is overlooked by strictly algebraic procedures in order to retrieve part of the loss efficiency and furthermore, to obtain insights that make them adaptable to different reservoir situations. In this work, we show that this physical information may be possible by aggregation of system coefficients via percolation. The process can be efficiently encapsulated into a two-stage preconditioner that relies on the solution of the identified highly connected blocks (i.e., aggregates) followed by a deflation preconditioner. The proposed approach is coined as two-stage percolation aggregation (2SPA) and emerges as a representative member of a new family of physics-driven solvers to exploit connectivity (and consequently, flow trends) on highly heterogeneous media. We compare the performance of 2SPA preconditioner against ILU and conclude that the method could be promising to tackle a wider range of reservoir scenarios.