  We consider the P1/P1 or P1b/P1 finite element approximations to the Stokes
equations in a bounded smooth domain subject to the slip boundary condition. A
penalty method is applied to address the essential boundary condition $u\cdot n
= g$ on $\partial\Omega$, which avoids a variational crime and simultaneously
facilitates the numerical implementation. We give $O(h^{1/2} + \epsilon^{1/2} +
h/\epsilon^{1/2})$-error estimate for velocity and pressure in the energy norm,
where $h$ and $\epsilon$ denote the discretization parameter and the penalty
parameter, respectively. In the two-dimensional case, it is improved to $O(h +
\epsilon^{1/2} + h^2/\epsilon^{1/2})$ by applying reduced-order numerical
integration to the penalty term. The theoretical results are confirmed by
numerical experiments.
