Summary Gradient reproducing kernel particle method (GRKPM) is a meshless technique which incorporates the first gradients of the function into the reproducing equation of RKPM. Therefore, in two-dimensionalspace GRKPM introducesthree types of shape functions rather than one. The robustness of GRKPM’s shape functions is established by reconstruction of a third-order polynomial. To enforce the essential boundary conditions (EBCs), GRKPM’s shape functions are modified by transformation technique. By utilizing the modified shape functions, the weak form of the nonlinear evolutionary Buckley-Leverett (BL) equation is discretized in space, rendering a system of nonlinear ordinary differential equations (ODEs). Subsequently, Gear’s method is applied for temporal discretization of the ODEs. Through numerical experiments, employment of a moderate viscosity seeks the efficacy ofthesolutionwhen thediffusionterm isimportant;moreover, applicationof a small viscosity confirms the potential of the approach for treatment of the problems involving steep gradient regions. The outcomes are verified by performing convergence tests using uniformly spaced particles. Consideration of non-uniform distribution of particles further demonstrates the virtue of the presented methodology in producing smooth profiles in the critical regions near the fronts.