The construction of valid and flexible cross-covariance functions is a fundamental task for modeling multivariate space-time data arising from climatological and oceanographical phenomena. Indeed, a suitable specification of the covariance structure allows to capture both the space-time dependencies between the observations and the development of accurate predictions. For data observed over large portions of planet Earth it is necessary to take into account the curvature of the planet. Hence the need for random field models defined over spheres across time. In particular, the associated covariance function should depend on the geodesic distance, which is the most natural metric over the spherical surface. In this work, we propose a flexible parametric family of matrix-valued covariance functions, with both marginal and cross structure being of the Gneiting type. We additionally introduce a different multivariate Gneiting model based on the adaptation of the latent dimension approach to the spherical context. Finally, we assess the performance of our models through the study of a bivariate space-time data set of surface air temperatures and precipitations.