The treatment of higher order perturbations of branes is considered using a covariant variational approach. This covariant variational approach brings to the forefront the geometric structure of the underlying perturbation theory, as opposed to a more commonly used ``direct approach,'' that ignores the variational origin. In addition, it offers a clear calculational advantage with respect to so-called ``gauge fixed'' treatments that distinguish tangential and normal modes, as it emphasizes the symmetries of the geometric models that describe the brane dynamics. We restrict our attention to a brane action that depends at most on first derivatives of the embedding functions of the world volume spanned by the brane in its evolution. We consider first and second variations of the action that describes the brane dynamics. The first variation produces the equations of motion, as is well known. In the second variation we derive the Jacobi equations for these kind of models, and we emphasize the role of the Hessian matrix. This is extended to third order in variations, first in a flat and then in a curved spacetime background. Further, we specialize to the relevant case of the Dirac-Nambu-Goto action that describes extremal branes. The proper setting of a covariant variational approach allows to go in principle to geometric models that depend on higher derivatives of the embedding functions, and higher order perturbations, with the due complications involved, but with a solid framework in place.