Owing to relatively rapid water content variations with respect to electromagnetic wavelength at GPR frequencies, the vadose zone usually exhibits significant velocity anisotropy. Neglecting anisotropy in traveltime tomographic reconstruction leads to artifacts that can obscure important subsurface features, the vadose zone being subject to this problem. In this paper, an algorithm for crosshole GPR geostatistical traveltime tomography in elliptically anisotropic media is presented. The advantages of the geostatistical tomography algorithm are that the solution is regularized by the covariance of the model parameters and that stochastic simulations can be performed to appraise the variability of the solution space. The implementation relies on a fast curved raytracing scheme specifically crafted for the problem at hand. The benefits of the algorithm to image the vadose zone are illustrated through a synthetic case that is representative of typical studies in quaternary geological settings. The results show that considering anisotropy yields better fit to the data at high ray angles and reduces reconstruction artifacts.