In this article we define a surface finite element method (SFEM) for the numerical solution of parabolic partial differential equations on hypersurfaces $\Gamma$ in $\mathbb R^{n+1}$. The key idea is based on the approximation of $\Gamma$ by a polyhedral surface $\Gamma_h$ consisting of a union of simplices (triangles for $n=2$, intervals for $n=1$) with vertices on $\Gamma$. A finite element space of functions is then defined by taking the continuous functions on $\Gamma_h$ which are linear affine on each simplex of the polygonal surface. We use surface gradients to define weak forms of elliptic operators and naturally generate weak formulations of elliptic and parabolic equations on $\Gamma$. Our finite element method is applied to weak forms of the equations. The computation of the mass and element stiffness matrices is simple and straightforward. We give an example of error bounds in the case of semi-discretization in space for a fourth order linear problem. Numerical experiments are described for several linear and nonlinear partial differential equations. In particular, the power of the method is demonstrated by employing it to solve highly nonlinear second and fourth order problems such as surface Allen-Cahn and Cahn-Hilliard equations and surface level set equations for geodesic mean curvature flow.