We consider the numerical approximation of high order Partial Differential Equations (PDEs) defined on surfaces in the three dimensional space, with particular emphasis on closed surfaces. We consider computational domains that can be represented by B-splines or NURBS, as for example the sphere, and we spatially discretize the PDEs by means of NURBS-based Isogeometric Analysis in the framework of the standard Galerkin method. We numerically solve benchmark Laplace-Beltrami problems of the fourth and sixth order, as well as the corresponding eigenvalue problems, with the goal of analyzing the role of the continuity of the NURBS basis functions on closed surfaces. In this respect, we show that the use of globally high order continuous basis functions, as allowed by the construction of periodic NURBS, leads to the efficient solution of the high order PDEs. Finally, we consider the numerical solution of high order phase field problems on closed surfaces, namely the Cahn-Hilliard and crystal equations.