Abstract: A proof of convergence is given for semi- and full discretizations of mean curvature flow of closed two-dimensional surfaces. The numerical method proposed and studied here combines evolving finite elements, whose nodes determine the discrete surface like in Dziuk’s method, and linearly implicit backward difference formulae for time integration. The proposed method differs from Dziuk’s approach in that it discretizes Huisken’s evolution equations for the normal vector and mean curvature and uses these evolving geometric quantities in the velocity law projected to the finite element space. This numerical method admits a convergence analysis in the case of finite elements of polynomial degree at least two and backward difference formulae up to order five. The error analysis combines stability estimates and consistency estimates to yield optimal-order $H^1$-norm error bounds for the computed surface position, velocity, normal vector and mean curvature. Numerical experiments illustrate and complement the theoretical results.