A new set of orthonormal vector polynomials in a square area, which can be used in image distortion mapping and wavefront gradient vector datum fitting, is derived. These vector polynomials are developed from the gradients of the circular Zernike polynomials orthonormalization by using Gram-Schmidt technique. When the slope is fitted by these vector polynomials, the fitting coefficients can be derived and transformed to the wavefront description of the Zernike polynomials mode by using a linear transform, and the phase information is then extracted. Experimental results show that the slope data from Shack-Hartmann wavefront sensor over a square area are well fitted by theses vector polynomials. The vector polynomial wavefront reconstruction method can reconstruct the tested wavefront quite well and achieve the same accuracy as Southwell zonal method does.