Multivariate spatial modeling is a rapidly growing field, but most extant models are infeasible for use with massive spatial processes. In this work we introduced a highly flexible, interpretable and scalable multiresolution approach to multivariate spatial modeling. Relying on compactly supported basis functions and Gaussian Markov random field specifications for coefficients results in efficient and scalable calculation routines for likelihood evaluations and cokriging. We analytically show that special parameterizations approximate popular existing models. Moreover, the multiresolution approach allows for arbitrary specification of scale dependence between processes. We illustrate our approach through Monte Carlo studies to illustrate implied stochastic behavior and test our ability to recover scale dependence, and moreover examine a complex large bivariate observational minimum and maximum temperature dataset over the western United States.