This paper deals with vector (or multivariate) random fields in space and/or time with second-order moments, for which a framework is needed for specifying not only the properties of each component but also the possible cross relationships among the components. We derive basic properties of the covariance matrix function of the vector random field and propose three approaches to construct covariance matrix functions for Gaussian or non-Gaussian random fields. The first approach is to take derivatives of a univariate covariance function, the second one is to work on the univariate random field whose index domain is in a higher dimension and the third one is based on the scale mixture of separable spatio-temporal covariance matrix functions. To illustrate these methods, many parametric or semiparametric examples are formulated.