  Matrix-valued covariance functions are crucial to geostatistical modeling of
multivariate spatial data. The classical assumption of symmetry of a
multivariate covariance function is overlay restrictive and has been considered
as unrealistic for most of real data applications. Despite of that, the
literature on asymmetric covariance functions has been very sparse. In
particular, there is some work related to asymmetric covariances on Euclidean
spaces, depending on the Euclidean distance. However, for data collected over
large portions of planet Earth, the most natural spatial domain is a sphere,
with the corresponding geodesic distance being the natural metric. In this
work, we propose a strategy based on spatial rotations to generate asymmetric
covariances for multivariate random fields on the $d$-dimensional unit sphere.
We illustrate through simulations as well as real data analysis that our
proposal allows to achieve improvements in the predictive performance in
comparison to the symmetric counterpart.
