The Laplace-Beltrami system of nonlinear, elliptic, partial differential equations has utility in the generation of computational grids on complex and highly curved geometry. Discretization of this system with Finite Elements readily accommodates unstructured grids, but generates a large, sparse, ill-conditioned system of nonlinear discrete equations. The extensive use of the Laplace-Beltrami approach, particularly in large-scale applications, has been limited by the scalability and efficiency of solvers. This paper addresses this limitation by developing two nonlinear solvers based on the Jacobian-Free Newton-Krylov (JFNK) methodology. A key feature of these methods is that the Jacobian is not formed explicitly for use by the underlying linear solver. Iterative linear solvers such as the Generalized Minimal RESidual (GMRES) method do not technically require the stand-alone Jacobian; instead its action on a vector is approximated through two nonlinear function evaluations. The preconditioning required by GMRES is also discussed; two different preconditioners are developed, both of which are readily treated with existing Algebraic Multigrid (AMG) methods. Further, the most efficient preconditioner overall for the problems considered is based on a Picard linearization. Numerical examples demonstrate that these new solvers are significantly faster than a standard Newton-Krylov approach; a speedup factor of approximately 26 was obtained for the Picard preconditioner on the largest grids studied here. In addition, these JFNK solvers exhibit good algorithmic scaling with increasing grid size.