This paper proposes a new class of covariance functions for bivariate random fields on spheres, having the same properties as the bivariate Matérn model proposed in Euclidean spaces. The new class depends on the geodesic distance on a sphere; it allows for indexing differentiability (in the mean square sense) and fractal dimensions of the components of any bivariate Gaussian random field having such covariance structure. We find parameter conditions ensuring positive definiteness. We discuss other possible models and illustrate our findings through a simulation study, where we explore the performance of maximum likelihood estimation method for the parameters of the new covariance function. A data illustration then follows, through a bivariate data set of temperatures and precipitations, observed over a large portion of the Earth, provided by the National Oceanic and Atmospheric Administration Earth System Research Laboratory.