This paper considers a multivariate spatial random field, with each component having univariate marginal distributions of the skew-Gaussian type. We assume that the field is defined spatially on the unit sphere embedded in R 3 , allowing for modeling data available over large portions of planet Earth. This model admits explicit expressions for the marginal and cross covariances. However, the n-dimensional distributions of the field are difficult to evaluate, because it requires the sum of 2 n terms involving the cumulative and probability density functions of a n-dimensional Gaussian distribution. Since in this case inference based on the full likelihood is computationally unfeasible, we propose a composite likelihood approach based on pairs of spatial observations. This last being possible thanks to the fact that we have a closed form expression for the bivariate distribution. We illustrate the effectiveness of the method through simulation experiments and the analysis of a real data set of minimum and maximum surface air temperatures.