We develop a finite element method for the vector Laplacian based on the covariant derivative of tangential vector fields on surfaces embedded in R3. Closely related operators arise in models of flow on surfaces as well as elastic membranes and shells. The method is based on standard continuous parametric Lagrange elements that describe a R3 vector field on the surface, and the tangent condition is weakly enforced using a penalization term. We derive error estimates that take into account the approximation of both the geometry of the surface and the solution to the partial differential equation. In particular, we note that to achieve optimal order error estimates, in both energy and L2 norms, the normal approximation used in the penalization term must be of the same order as the approximation of the solution. This can be fulfilled either by using an improved normal in the penalization term, or by increasing the order of the geometry approximation. We also present numerical results using higher-order finite elements that verify our theoretical findings.