A novel phase field model for Willmore flow is proposed based on a nested variational time discretization. Thereby, the mean curvature in the Willmore functional is replaced by an approximate speed of mean curvature motion, which is computed via a fully implicit variational model for time discrete mean curvature motion. The time discretization of Willmore flow is then performed in a nested fashion: in an outer variational approach a natural time discretization is setup for the actual Willmore flow, whereas for the involved mean curvature the above variational approximation is taken into account. Hence, in each time step a PDE-constrained optimization problem has to be solved in which the actual surface geometry as well as the geometry resulting from the implicit curvature motion time step are represented by phase field functions. The convergence behavior is experimentally validated and compared with rigorously proved convergence estimates for a simple linear model problem. Computational results in 2D and 3D underline the robustness of the new discretization, in particular for large time steps and in comparison with a semi-implicit convexity splitting scheme. Furthermore, the new model is applied as a minimization method for elastic functionals in image restoration.