A new numerical technique based on radial basis functions (RBFs) is presented for fitting a vector field tangent to the sphere, S 2 , from samples of the field at "scattered" locations on S 2 . The method naturally provides a way to decompose the reconstructed field into its individual Helmholtz-Hodge components, i.e. into divergence-free and curl-free parts, which is useful in many applications from the atmospheric and oceanic sciences (e.g., in diagnosing the horizontal wind and ocean currents). Several approximation results for the method will be derived. In particular, Sobolev-type error estimates are obtained for both the interpolant and its decomposition. Optimal stability estimates for the associated interpolation matrices are also presented. Finally, numerical validation of the theoretical results is given for vector fields with similar characteristics to those of atmospheric wind fields.