A new solving method for Laplace equation with over-determined geodetic boundary conditions is proposed in the paper, with the help of minimizing some kinds of quadratic functional in calculus of variation. At first, the so-called variational solution for over-determined geodetic boundary value problem is defined in terms of principles in calculus of variation. Then theoretical properties related with the solution are derived, especially for its existence, uniqueness and optimal approximation. And then the computational method of the solution is discussed, and its expression is exhibited under the case that all boundaries are spheres. Finally an arithmetic example about EGM96 gravity field model is given, and the computational resulta show that the proposed method can efficiently raise accuracy to deal with gravity data. In all, the variational solution of over-determined geodetic boundary value problem can not only fit to deal with many kinds of gravity data in a united form, but also has strict mathematical basements.