One major challenge in modeling Darcy flow in heterogeneous porous media is simulating the material interfaces accurately. To overcome this defect, the refraction law is fully introduced into the numerical manifold method (NMM) as an a posteriori condition. To achieve a better accuracy of the Darcy velocity and continuous nodal velocity, a high-order weight function with a continuous nodal gradient is adopted. NMM is an advanced method with two independent cover systems, which can easily solve both continuous and discontinuous problems in a unified form. Moreover, a regular mathematical mesh, independent of the physical domain, is used in the NMM model. Compared to the conforming mesh of other numerical methods, it is more efficient and flexible. A number of numerical examples were simulated by the new NMM model, comparing the results with the original NMM model and the analytical solutions. Thereby, it is proven that the proposed method is accurate, efficient, and robust for modeling Darcy flow in heterogeneous porous media, while the refraction law is satisfied rigorously.