Modeling the variability of brain structures is a fundamental problem in the neurosciences. In this paper, we start from a dataset of precisely delineated anatomical structures in the cerebral cortex: a set of 72 sulcal lines in each of 98 healthy human subjects. We propose an original method to compute the average sulcal curves, which constitute the mean anatomy in this context. The second order moment of the sulcal distribution is modeled as a sparse field of covariance tensors (symmetric, positive definite matrices). To extrapolate this information to the full brain, one has to overcome the limitations of the standard Euclidean matrix calculus. We propose an affine-invariant Riemannian framework to perform computations with tensors. In particular, we generalize radial basis function (RBF) interpolation and harmonic diffusion PDEs to tensor fields. As a result, we obtain a dense 3D variability map which proves to be in accordance with previously published results on smaller samples s...