We study numerical approximations for geometric evolution equations arising as gradient flows for energy functionals that are quadratic in the principal curvatures of a two-dimensional surface. Beside the well-known Willmore and Helfrich flows we will also consider flows involving the Gaussian curvature of the surface. Boundary conditions for these flows are highly nonlinear, and we use a variational approach to derive weak formulations, which naturally can be discretized with the help of a mixed finite element method. Our approach uses a parametric finite element method,
which can be shown to lead to good mesh properties. We prove stability estimates for a semidiscrete (discrete in space, continuous in time) version of the method and show
existence and uniqueness results in the fully discrete case. Finally, several numerical results are presented involving convergence tests as well as the first computations with Gaussian curvature and/or free or semi-free boundary conditions.