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 subjects. Moreover, leave one (sulcus) out tests show that our model is globally able to recover the missing information when there is a consistent neighboring variability. Last but not least, we propose innovative methods to analyze the asymmetry of brain variability. As expected, the greatest asymmetries are found in regions that includes the primary language areas. Interestingly, such an asymmetry in anatomical variance could explain why there may be greater power to detect group activation in one hemisphere than the other in fMRI studies.