Field trials for variety selection often exhibit spatial correlation between plots. When multivariate data are analysed from these field trials, there is the added complication in having to simultaneously account for correlation between the traits at both the residual and genetic levels. This may be temporal correlation in the case of multi-harvest data from perennial crop field trials, or between-trait correlation in multi-trait data sets. Use of parsimonious yet plausible models for the variance-covariance structure of the residuals for such data is a key element to achieving an efficient and inferentially sound analysis. In this paper, a model is developed for the residual variance-covariance structure firstly by considering a multivariate autoregressive model in one spatial direction and then extending this to two spatial directions. Conditions for ensuring that the processes are directionally invariant are presented. Using a canonical decomposition, these directionally invariant processes can be transformed into a set of independent separable processes. This simplifies the estimation process. The new model allows for flexible modelling of the spatial and multivariate interaction and allows for different spatial correlation parameters for each harvest or trait. The methods are illustrated using data from lucerne breeding trials at several environments.