In this paper, the variational multiscale interpolating element-free Galerkin (VMIEFG) method is developed to obtain the numerical solution ofthenonlinearDarcy–Forchheimer model. We use the interpolating moving least squares method instead of the moving least squares approximation to construct meshless shape functions with delta function properties. Then the flux boundary condition of the Darcy–Forchheimer model can be handled easily. Hughes’ variational multiscale (HVM) method is applied to overcome the numerical oscillation caused by equal-order basis for the velocity and pressure. Moreover, the HVM ensures that the resultant formulation in the VMIEFG method is consistent and the stabilization parameter (or tensor) appears naturally. Consequently, the stabilization parameter is free of user-defined. The fixed point iteration method is used to deal with the nonlinear term. Some numerical examples are provided to illustrate the stability and performance of the proposed method for solving the nonlinear Darcy–Forchheimer model.