In this paper we first revisit a classical problem of computing variational splines. We propose to compute local variational splines in the sense that they are interpolatory splines which minimize the energy norm over a subinterval. We shall show that the error between local and global variational spline interpolants decays exponentially over a fixed subinterval as the support of the local variational spline increases. By piecing together these locally defined splines, one can obtain a very good C0 approximation of the global variational spline. Finally we generalize this idea to approximate global tensor product B-spline interpolatory surfaces.