Multivariate space-time data are increasingly available in various scientific disciplines. When analyzing these data, one of the key issues is to describe the multivariate space-time dependencies. Under the Gaussian framework, one needs to propose relevant models for multivariate space-time covariance functions, i.e. matrix-valued mappings with the additional requirement of non-negative definiteness. We propose a flexible parametric class of cross-covariance functions for multivariate space-time Gaussian random fields. Space-time components belong to the (univariate) Gneiting class of space-time covariance functions, with Matern or Cauchy covariance functions in the spatial margins. The smoothness and scale parameters can be different for each variable. We provide sufficient conditions for positive definiteness. A simulation study shows that the parameters of this model can be efficiently estimated using weighted pairwise likelihood, which belongs to the class of composite likelihood methods. We then illustrate the model on a French dataset of weather variables.