In the recent years, there has been a growing interest in proposing covariance models for multivariate Gaussian random fields. Some of these covariance models are very flexible and can capture both the marginal and the cross-spatial dependence of the components of the associated multivariate Gaussian random field. However, effective estimation methods for these models are somehow unexplored. Maximum likelihood is certainly a useful tool, but it is impractical in all the circumstances where the number of observations is very large. In this work, we consider two possible approaches based on composite likelihood for multivariate covariance model estimation. We illustrate, through simulation experiments, that our methods offer a good balance between statistical efficiency and computational complexity. Asymptotic properties of the proposed estimators are assessed under increasing domain asymptotics. Finally, we apply the method for the analysis of a bivariate dataset on chlorophyll concentration and sea surface temperature in the Chilean coast. © 2016, International Biometric Society.