We consider the penalty method for the stationary Navier-Stokes equations with the slip boundary condition. The well-posedness and the regularity theorem of the penalty problem are investigated, and we obtain the optimal error estimate O( ) in H k -norm, where is the penalty parameter. We are concerned with the finite element approximation with the P1b/P1 element to the penalty problem. The well-posedness of discrete problem is proved. We obtain the error estimate O(h + √ + h/ √ ) for the non-reduced-integration scheme with d = 2, 3, and the reduced-integration scheme with d = 3, where h is the discretization parameter and d is the spatial dimension. For the reduced-integration scheme with d = 2, we prove the convergence order O(h + √ + h 2 / √ ). The theoretical results are verified by numerical experiments.