Joint modeling of spatially-oriented dependent variables are commonplace in the environmental sciences, where scientists seek to estimate the relationships among a set of environmental outcomes accounting for dependence among these outcomes and the spatial dependence for each outcome. Such modeling is now sought for very large data sets where the variables have been measured at a very large number of locations. Bayesian inference, while attractive for accommodating uncertainties through their hierarchical structures, can become computationally onerous for modeling massive spatial data sets because of their reliance on iterative estimation algorithms. This manuscript develops a conjugate Bayesian framework for analyzing multivariate spatial data using analytically tractable posterior distributions that do not require iterative algorithms. We discuss differences between modeling the multivariate response itself as a spatial process and that of modeling a latent process. We illustrate the computational and inferential benefits of these models using simulation studies and real data analyses for a Vege Indices dataset with observed locations numbering in the millions.