This paper deals with vector field interpolation, i.e., the data are R-3 values located in scattered R-3 points, while the interpolating function is a function from R-3 into R-3. In order to take into account possible connections between the components of the interpolant, we derive it by solving a variational spline problem involving the rotational and the divergence of the interpolant, and depending on a parameter rho significative of the balance of the rotational part and of the divergence part, and on the order m of derivatives of the rotational and divergence involved in the minimized seminorm. We so obtain interpolants whose expression is sigma(x) = Sigma(i=1)(n) Phi(x - x(i))a(i) + p(m-1) (x), where Phi is some 3 x 3 matrix function, p(m-1) is a degree m-1 vectorial polynomial, and where the a(i) are R-3-vectors. Besides, the a(i) meet a relation generalizing the usual orthogonality to all polynomials of degree at most m-1. For rho = 1, we find the usual m-harmonic splines in each component of sigma. Numerical examples show the interest of the method, and we compare the so-obtained functions with the ones obtained by Matlab's procedures. (C) 2002 Elsevier Science Ltd. All rights reserved.