We describe our implementation of the multivariate Matern model for multivariate spatial datasets, using Vecchia's approximation and a Fisher scoring optimization algorithm. We consider various pararameterizations for the multivariate Mat\'ern that have been proposed in the literature for ensuring model validity, as well as an unconstrained model. A strength of our study is that the code is tested on many real-world multivariate spatial datasets. We use it to study the effect of ordering and conditioning in Vecchia's approximation and the restrictions imposed by the various parameterizations. We also consider a model with nuggets that are correlated across components and find that forcing the nugget correlation to be zero can have a serious impact on the other model parameters, so we suggest allowing correlated nuggets.