Abstract: This paper proposes a gradient–divergence (grad–div) stabilized two-level finite element method. The method adopts a hierarchical strategy: first, a nonlinear Navier-Stokes-type variational inequality problem incorporating a grad–div stabilization term is solved on a coarse grid to yield a coarse-grid solution; subsequently, the nonlinear convective term is linearized via the coarse-grid solution on a fine grid, and a linearized problem (of Oseen, Stokes, or Newton type) is solved to derive the final solution. The grad–div stabilization term serves to suppress pressure's influence on velocity. Theoretical error bounds for the approximate solutions of this numerical method are established, and numerical experiments validate the correctness of the theoretical results as well as the method's effectiveness. Numerical findings demonstrate that, in comparison with the standard two-level finite element method, the proposed method achieves a significant enhancement in the accuracy of velocity approximations.